source: branches/stable/GDE/MAFFT/mafft-7.055-with-extensions/core/io.c

Last change on this file was 10371, checked in by aboeckma, 11 years ago

updated mafft version. Added extensions (no svn ignore, yet)

File size: 115.1 KB
Line 
1#include "mltaln.h"
2
3static int upperCase = 0;
4
5#define DEBUG   0
6#define IODEBUG 0
7
8char 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
55void 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
64void 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
78static 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
212static 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}
336void 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}
481void 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
612void 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
661void 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}
790void 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
917char *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
935void ErrorExit( char *message )
936{
937        fprintf( stderr, "%s\n", message );
938        exit( 1 );
939}
940
941void 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       
951void seqUpper( int nseq, char **seq ) /* not used */
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
962void 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
973int 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
983int getaline_fp_eof_new(s, l, fp)  /* end of file -> return 1 */
984char    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
999int myfgets(s, l, fp)  /* l°ÊŸå€Ï¡¢¹ÔËö€Þ€ÇÆÉ€ßÈô€Ð€¹ */
1000char    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
1015float 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
1032void 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
1063int allSpace( char *str )
1064{
1065        int value = 1;
1066        while( *str ) value *= ( !isdigit( *str++ ) );
1067        return( value );
1068}
1069       
1070void 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
1077void 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
1186static 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
1203void 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
1212static 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 == '>' )
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, "=== '>' in sequence is unacceptable.\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
1240static int onlyAlpha_lower( char *str )
1241{
1242        char tmp;
1243        char *res = str;
1244        char *bk = str;
1245
1246        while( (tmp=*str++) )
1247                if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1248                        *res++ = tolower( tmp );
1249        *res = 0;
1250        return( res - bk );
1251}
1252static int onlyAlpha_upper( char *str )
1253{
1254        char tmp;
1255        char *res = str;
1256        char *bk = str;
1257
1258        while( (tmp=*str++) )
1259                if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1260                        *res++ = toupper( tmp );
1261        *res = 0;
1262        return( res - bk );
1263}
1264
1265void kake2hiku( char *str )
1266{
1267        do
1268                if( *str == '*' ) *str = '-';
1269        while( *str++ );
1270}
1271
1272char *load1SeqWithoutName_realloc_casepreserve( FILE *fpp )
1273{
1274        int c, b;
1275        char *cbuf;
1276        int size = N;
1277        char *val;
1278
1279        val = malloc( (size+1) * sizeof( char ) );
1280        cbuf = val;
1281
1282        b = '\n';
1283        while( ( c = getc( fpp ) ) != EOF &&           
1284          !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1285        {
1286                *cbuf++ = (char)c;  /* Ĺ€¹€®€Æ€â€·€é€Ê€€ */
1287                if( cbuf - val == size )
1288                {
1289                        size += N;
1290                        fprintf( stderr, "reallocating...\n" );
1291                        val = (char *)realloc( val, (size+1) * sizeof( char ) );
1292                        if( !val )
1293                        {
1294                                fprintf( stderr, "Allocation error in load1SeqWithoutName_realloc \n" );
1295                                exit( 1 );
1296                        }
1297                        fprintf( stderr, "done.\n" );
1298                        cbuf = val + size-N;
1299                }
1300                b = c;
1301        }
1302        ungetc( c, fpp );
1303        *cbuf = 0;
1304        onlyGraph( val );
1305//      kake2hiku( val );
1306        return( val );
1307}
1308
1309char *load1SeqWithoutName_realloc( FILE *fpp )
1310{
1311        int c, b;
1312        char *cbuf;
1313        int size = N;
1314        char *val;
1315
1316        val = malloc( (size+1) * sizeof( char ) );
1317        cbuf = val;
1318
1319        b = '\n';
1320        while( ( c = getc( fpp ) ) != EOF &&           
1321          !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1322        {
1323                *cbuf++ = (char)c;  /* Ĺ€¹€®€Æ€â€·€é€Ê€€ */
1324                if( cbuf - val == size )
1325                {
1326                        size += N;
1327                        fprintf( stderr, "reallocating...\n" );
1328                        val = (char *)realloc( val, (size+1) * sizeof( char ) );
1329                        if( !val )
1330                        {
1331                                fprintf( stderr, "Allocation error in load1SeqWithoutName_realloc \n" );
1332                                exit( 1 );
1333                        }
1334                        fprintf( stderr, "done.\n" );
1335                        cbuf = val + size-N;
1336                }
1337                b = c;
1338        }
1339        ungetc( c, fpp );
1340        *cbuf = 0;
1341        if( dorp == 'd' )
1342                onlyAlpha_lower( val );
1343        else
1344                onlyAlpha_upper( val );
1345        kake2hiku( val );
1346        return( val );
1347}
1348
1349int load1SeqWithoutName_new( FILE *fpp, char *cbuf )
1350{
1351        int c, b;
1352        char *bk = cbuf;
1353
1354        b = '\n';
1355        while( ( c = getc( fpp ) ) != EOF &&                    /* by T. Nishiyama */
1356          !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1357        {
1358                *cbuf++ = (char)c;  /* Ĺ€¹€®€Æ€â€·€é€Ê€€ */
1359                b = c;
1360        }
1361        ungetc( c, fpp );
1362        *cbuf = 0;
1363        if( dorp == 'd' )
1364                onlyAlpha_lower( bk );
1365        else
1366                onlyAlpha_upper( bk );
1367        kake2hiku( bk );
1368        return( 0 );
1369}
1370
1371
1372void readDataforgaln( FILE *fp, char **name, int *nlen, char **seq )
1373{
1374        int i; 
1375        static char *tmpseq = NULL;
1376
1377#if 0
1378        if( !tmpseq )
1379        {
1380                tmpseq = AllocateCharVec( N );
1381        }
1382#endif
1383
1384        rewind( fp );
1385        searchKUorWA( fp );
1386
1387        for( i=0; i<njob; i++ )
1388        {
1389                name[i][0] = '='; getc( fp ); 
1390#if 0
1391                fgets( name[i]+1, B-2, fp );
1392                j = strlen( name[i] );
1393                if( name[i][j-1] != '\n' )
1394                        ErrorExit( "Too long name\n" );
1395                name[i][j-1] = 0;
1396#else
1397                myfgets( name[i]+1, B-2, fp ); 
1398#endif
1399#if 0
1400                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1401#endif
1402                tmpseq = load1SeqWithoutName_realloc( fp );
1403                strcpy( seq[i], tmpseq );
1404                nlen[i] = strlen( seq[i] );
1405                free( tmpseq );
1406        }
1407        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1408#if 0
1409        free( tmpseq );
1410#endif
1411}
1412
1413void readData_varlen( FILE *fp, char **name, int *nlen, char **seq )
1414{
1415        int i; 
1416        static char *tmpseq = NULL;
1417
1418        rewind( fp );
1419        searchKUorWA( fp );
1420
1421        for( i=0; i<njob; i++ )
1422        {
1423                name[i][0] = '='; getc( fp ); 
1424#if 0
1425                fgets( name[i]+1, B-2, fp );
1426                j = strlen( name[i] );
1427                if( name[i][j-1] != '\n' )
1428                        ErrorExit( "Too long name\n" );
1429                name[i][j-1] = 0;
1430#else
1431                myfgets( name[i]+1, B-2, fp ); 
1432#endif
1433#if 0
1434                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1435#endif
1436                tmpseq = load1SeqWithoutName_realloc( fp );
1437                nlen[i] = strlen( tmpseq );
1438//              fprintf( stderr, "nlen[%d] = %d\n", i+1, nlen[i] );
1439                seq[i] = calloc( nlen[i]+1, sizeof( char ) );
1440                strcpy( seq[i], tmpseq );
1441                free( tmpseq );
1442        }
1443        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1444#if 0
1445        free( tmpseq );
1446#endif
1447}
1448
1449void readData_pointer2( FILE *fp, int nseq, char **name, int *nlen, char **seq )
1450{
1451        int i; 
1452        static char *tmpseq = NULL;
1453
1454#if 0
1455        if( !tmpseq )
1456        {
1457                tmpseq = AllocateCharVec( N );
1458        }
1459#endif
1460
1461        rewind( fp );
1462        searchKUorWA( fp );
1463
1464        for( i=0; i<nseq; i++ )
1465        {
1466                name[i][0] = '='; getc( fp ); 
1467#if 0
1468                fgets( name[i]+1, B-2, fp );
1469                j = strlen( name[i] );
1470                if( name[i][j-1] != '\n' )
1471                        ErrorExit( "Too long name\n" );
1472                name[i][j-1] = 0;
1473#else
1474                myfgets( name[i]+1, B-2, fp ); 
1475#endif
1476#if 0
1477                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1478#endif
1479                tmpseq = load1SeqWithoutName_realloc( fp );
1480                strcpy( seq[i], tmpseq );
1481                free( tmpseq );
1482                nlen[i] = strlen( seq[i] );
1483        }
1484        if( dorp == 'd' && upperCase != -1 ) seqLower( nseq, seq );
1485#if 0
1486        free( tmpseq );
1487#endif
1488        if( outnumber )
1489        {
1490                char *namebuf;
1491                char *cptr;
1492                namebuf = calloc( B+100, sizeof( char ) );
1493                for( i=0; i<nseq; i++ )
1494                {
1495                        namebuf[0] = '=';
1496                        cptr = strstr( name[i], "_numo_e_" );
1497                        if( cptr )
1498                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, cptr+8 );
1499                        else
1500                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, name[i]+1 );
1501                        strncpy( name[i], namebuf, B );
1502                        name[i][B-1] = 0;
1503                }
1504                free( namebuf );
1505//              exit( 1 );
1506        }
1507}
1508
1509
1510void readData_pointer_casepreserve( FILE *fp, char **name, int *nlen, char **seq )
1511{
1512        int i; 
1513        static char *tmpseq = NULL;
1514
1515#if 0
1516        if( !tmpseq )
1517        {
1518                tmpseq = AllocateCharVec( N );
1519        }
1520#endif
1521
1522        rewind( fp );
1523        searchKUorWA( fp );
1524
1525        for( i=0; i<njob; i++ )
1526        {
1527                name[i][0] = '='; getc( fp ); 
1528#if 0
1529                fgets( name[i]+1, B-2, fp );
1530                j = strlen( name[i] );
1531                if( name[i][j-1] != '\n' )
1532                        ErrorExit( "Too long name\n" );
1533                name[i][j-1] = 0;
1534#else
1535                myfgets( name[i]+1, B-2, fp ); 
1536#endif
1537#if 0
1538                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1539#endif
1540                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1541                strcpy( seq[i], tmpseq );
1542                free( tmpseq );
1543                nlen[i] = strlen( seq[i] );
1544        }
1545}
1546
1547
1548void readData_pointer( FILE *fp, char **name, int *nlen, char **seq )
1549{
1550        int i; 
1551        static char *tmpseq = NULL;
1552
1553#if 0
1554        if( !tmpseq )
1555        {
1556                tmpseq = AllocateCharVec( N );
1557        }
1558#endif
1559
1560        rewind( fp );
1561        searchKUorWA( fp );
1562
1563        for( i=0; i<njob; i++ )
1564        {
1565                name[i][0] = '='; getc( fp ); 
1566#if 0
1567                fgets( name[i]+1, B-2, fp );
1568                j = strlen( name[i] );
1569                if( name[i][j-1] != '\n' )
1570                        ErrorExit( "Too long name\n" );
1571                name[i][j-1] = 0;
1572#else
1573                myfgets( name[i]+1, B-2, fp ); 
1574#endif
1575#if 0
1576                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1577#endif
1578                tmpseq = load1SeqWithoutName_realloc( fp );
1579                strcpy( seq[i], tmpseq );
1580                free( tmpseq );
1581                nlen[i] = strlen( seq[i] );
1582        }
1583        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1584#if 0
1585        free( tmpseq );
1586#endif
1587        if( outnumber )
1588        {
1589                char *namebuf;
1590                char *cptr;
1591                namebuf = calloc( B+100, sizeof( char ) );
1592                for( i=0; i<njob; i++ )
1593                {
1594                        namebuf[0] = '=';
1595                        cptr = strstr( name[i], "_numo_e_" );
1596                        if( cptr )
1597                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, cptr+8 );
1598                        else
1599                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, name[i]+1 );
1600                        strncpy( name[i], namebuf, B );
1601                        name[i][B-1] = 0;
1602                }
1603                free( namebuf );
1604//              exit( 1 );
1605        }
1606}
1607
1608void readData( FILE *fp, char name[][B], int nlen[], char **seq )
1609{
1610        int i; 
1611        static char *tmpseq = NULL;
1612
1613#if 0
1614        if( !tmpseq )
1615        {
1616                tmpseq = AllocateCharVec( N );
1617        }
1618#endif
1619
1620        rewind( fp );
1621        searchKUorWA( fp );
1622
1623        for( i=0; i<njob; i++ )
1624        {
1625                name[i][0] = '='; getc( fp ); 
1626#if 0
1627                fgets( name[i]+1, B-2, fp );
1628                j = strlen( name[i] );
1629                if( name[i][j-1] != '\n' )
1630                        ErrorExit( "Too long name\n" );
1631                name[i][j-1] = 0;
1632#else
1633                myfgets( name[i]+1, B-2, fp ); 
1634#endif
1635#if 0
1636                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1637#endif
1638                tmpseq = load1SeqWithoutName_realloc( fp );
1639                strcpy( seq[i], tmpseq );
1640                nlen[i] = strlen( seq[i] );
1641                free( tmpseq );
1642        }
1643        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1644#if 0
1645        free( tmpseq );
1646#endif
1647}
1648
1649void cutAlignment( FILE *fp, int **regtable, char **revtable, int *outtable, char **name, char **outseq )
1650{
1651        int i, j; 
1652        int outlen;
1653        static char *tmpseq = NULL;
1654        static char *dumname = NULL;
1655        char *fs, *rs;
1656        int npos, lpos;
1657        int startpos, endpos, seqlen;
1658
1659        if( dumname == NULL )
1660        {
1661                dumname = AllocateCharVec( N );
1662        }
1663
1664        rewind( fp );
1665        searchKUorWA( fp );
1666
1667
1668        npos = 0;
1669        for( i=0; i<njob; i++ )
1670        {
1671                dumname[0] = '>'; getc( fp ); 
1672                myfgets( dumname+1, B-1, fp ); 
1673                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1674
1675                if( outtable[i] )
1676                {
1677//                      putc( '>', stdout );
1678//                      puts( dumname+1 );
1679
1680       
1681                        strncat( name[npos], dumname, B-1 );
1682                        name[npos][B-1] = 0;
1683       
1684                        if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1685                        seqlen = strlen( tmpseq );
1686                        lpos = 0;
1687                        for( j=0; j<5; j++ )
1688                        {
1689                                if( regtable[0][j*2] == -1 && regtable[0][j*2+1] == -1 ) continue;
1690
1691                                startpos = regtable[0][j*2];
1692                                endpos   = regtable[0][j*2+1];
1693                                if( startpos > endpos )
1694                                {
1695                                        endpos   = regtable[0][j*2];
1696                                        startpos = regtable[0][j*2+1];
1697                                }
1698
1699                                if( startpos < 0 ) startpos = 0;
1700                                if( endpos   < 0 ) endpos   = 0;
1701                                if( endpos   >= seqlen ) endpos   = seqlen-1;
1702                                if( startpos >= seqlen ) startpos = seqlen-1;
1703
1704//                              fprintf( stderr, "startpos = %d, endpos = %d\n", startpos, endpos );
1705
1706                                outlen = endpos - startpos+1;
1707                                if( revtable[0][j] == 'f' )
1708                                {
1709//                                      fprintf( stderr, "regtable[%d][st] = %d\n", i, regtable[0][j*2+0] );
1710//                                      fprintf( stderr, "regtable[%d][en] = %d\n", i, regtable[0][j*2+1] );
1711//                                      fprintf( stderr, "outlen = %d\n", outlen );
1712//                                      fprintf( stdout, "%.*s\n", outlen, tmpseq+regtable[0][j*2] );
1713                                        strncpy( outseq[npos] + lpos, tmpseq+startpos, outlen );
1714                                        lpos += outlen;
1715                                }
1716                                else
1717                                {
1718                                        fs = AllocateCharVec( outlen+1 );
1719                                        rs = AllocateCharVec( outlen+1 );
1720
1721                                        fs[outlen] = 0;
1722                                        strncpy( fs, tmpseq+startpos, outlen );
1723                                        sreverse( rs, fs );
1724//                                      fprintf( stdout, "%s\n", rs );
1725                                        strncpy( outseq[npos] + lpos, rs, outlen );
1726                                        lpos += outlen;
1727                                        free( fs );
1728                                        free( rs );
1729                                }
1730                                outseq[npos][lpos] = 0;
1731                        }
1732                        npos++;
1733                }
1734                free( tmpseq );
1735        }
1736}
1737
1738void cutData( FILE *fp, int **regtable, char **revtable, int *outtable )
1739{
1740        int i, j; 
1741        int outlen, seqlen, startpos, endpos;
1742        static char *tmpseq = NULL;
1743        static char *dumname = NULL;
1744        char *fs, *rs;
1745
1746        if( dumname == NULL )
1747        {
1748                dumname = AllocateCharVec( N );
1749        }
1750
1751        rewind( fp );
1752        searchKUorWA( fp );
1753
1754        for( i=0; i<njob; i++ )
1755        {
1756                dumname[0] = '='; getc( fp ); 
1757                myfgets( dumname+1, B-2, fp ); 
1758                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1759
1760                if( outtable[i] )
1761                {
1762                        gappick_samestring( tmpseq );
1763                        putc( '>', stdout );
1764                        puts( dumname+1 );
1765
1766                        seqlen = strlen( tmpseq );
1767
1768                        if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1769                        if( outtable[i] == 2 )
1770                        {
1771                                startpos = 0;
1772                                endpos   = seqlen-1;
1773                                outlen = endpos - startpos + 1;
1774                                fprintf( stdout, "%.*s\n", outlen, tmpseq+startpos );
1775                        }
1776                        else
1777                        {
1778                                for( j=0; j<5; j++ )
1779                                {
1780                                        if( regtable[i][j*2] == -1 && regtable[i][j*2+1] == -1 ) continue;
1781       
1782                                        startpos = regtable[i][j*2];
1783                                        endpos   = regtable[i][j*2+1];
1784       
1785                                        if( startpos > endpos )
1786                                        {
1787                                                endpos   = regtable[i][j*2];
1788                                                startpos = regtable[i][j*2+1];
1789                                        }
1790       
1791                                        if( startpos < 0 ) startpos = 0;
1792                                        if( endpos   < 0 ) endpos   = 0;
1793                                        if( endpos   >= seqlen ) endpos   = seqlen-1;
1794                                        if( startpos >= seqlen ) startpos = seqlen-1;
1795       
1796                                        outlen = endpos - startpos + 1;
1797                                        if( revtable[i][j] == 'f' )
1798                                        {
1799                                                fprintf( stderr, "startpos = %d\n", startpos );
1800                                                fprintf( stderr, "endpos   = %d\n", endpos );
1801                                                fprintf( stderr, "outlen = %d\n", outlen );
1802                                                fprintf( stdout, "%.*s\n", outlen, tmpseq+startpos );
1803                                        }
1804                                        else
1805                                        {
1806                                                fs = AllocateCharVec( outlen+1 );
1807                                                rs = AllocateCharVec( outlen+1 );
1808       
1809                                                fs[outlen] = 0;
1810                                                strncpy( fs, tmpseq+startpos, outlen );
1811                                                sreverse( rs, fs );
1812                                                fprintf( stdout, "%s\n", rs );
1813                                                free( fs );
1814                                                free( rs );
1815                                        }
1816                                }
1817                        }
1818                }
1819                free( tmpseq );
1820        }
1821}
1822
1823void catData( FILE *fp )
1824{
1825        int i; 
1826        static char *tmpseq = NULL;
1827        static char *dumname = NULL;
1828//      char *cptr;
1829
1830        if( dumname == NULL )
1831        {
1832                dumname = AllocateCharVec( N );
1833        }
1834
1835        rewind( fp );
1836        searchKUorWA( fp );
1837
1838        for( i=0; i<njob; i++ )
1839        {
1840                dumname[0] = '='; getc( fp ); 
1841                myfgets( dumname+1, B-2, fp ); 
1842                if( outnumber )
1843                {
1844                        fprintf( stdout, ">_numo_s_%08d_numo_e_", i+1 );
1845                }
1846                else
1847                {
1848                        putc( '>', stdout );
1849                }
1850                puts( dumname+1 );
1851                tmpseq = load1SeqWithoutName_realloc( fp );
1852                if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1853                puts( tmpseq );
1854                free( tmpseq );
1855        }
1856}
1857
1858int countATGC( char *s, int *total )
1859{
1860        int nATGC;
1861        int nChar;
1862        char c;
1863        nATGC = nChar = 0;
1864
1865        if( *s == 0 ) 
1866        {
1867                *total = 0;
1868                return( 0 );
1869        }
1870
1871        do
1872        {
1873                c = tolower( *s );
1874                if( isalpha( c ) )
1875                {
1876                        nChar++;
1877                        if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1878                                nATGC++;
1879                }
1880        }
1881        while( *++s );
1882
1883        *total = nChar;
1884        return( nATGC );
1885}
1886
1887double countATGCbk( char *s )
1888{
1889        int nATGC;
1890        int nChar;
1891        char c;
1892        nATGC = nChar = 0;
1893
1894        do
1895        {
1896                c = tolower( *s );
1897                if( isalpha( c ) )
1898                {
1899                        nChar++;
1900                        if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1901                                nATGC++;
1902                }
1903        }
1904        while( *++s );
1905        return( (double)nATGC / nChar );
1906}
1907
1908
1909int countnogaplen( char *seq )
1910{
1911        int val = 0;
1912        while( *seq )
1913                if( *seq++ != '-' ) val++;
1914        return( val );
1915}
1916
1917int countnormalletters( char *seq, char *ref )
1918{
1919        int val = 0;
1920        while( *seq )
1921                if( strchr( ref, *seq++ ) ) val++;
1922        return( val );
1923}
1924
1925void getnumlen_casepreserve( FILE *fp, int *nlenminpt )
1926{
1927        int total;
1928        int nsite = 0;
1929        int atgcnum;
1930        int i, tmp;
1931        char *tmpseq, *tmpname;
1932        double atgcfreq;
1933        tmpname = AllocateCharVec( N );
1934        njob = countKUorWA( fp );
1935        searchKUorWA( fp );
1936        nlenmax = 0;
1937        *nlenminpt = 99999999;
1938        atgcnum = 0;
1939        total = 0;
1940        for( i=0; i<njob; i++ )
1941        {
1942                myfgets( tmpname, N-1, fp );
1943                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1944                tmp = strlen( tmpseq );
1945                if( tmp > nlenmax ) nlenmax  = tmp;
1946                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
1947                atgcnum += countATGC( tmpseq, &nsite );
1948                total += nsite;
1949                free( tmpseq );
1950        }
1951        free( tmpname );
1952        atgcfreq = (double)atgcnum / total;
1953//      fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1954        if( dorp == NOTSPECIFIED )
1955        {
1956                if( atgcfreq > 0.75 )   
1957                {
1958                        dorp = 'd';
1959                        upperCase = -1;
1960                }
1961                else                 
1962                {
1963                        dorp = 'p';
1964                        upperCase = 0;
1965                }
1966        }
1967}
1968
1969void getnumlen_nogap( FILE *fp, int *nlenminpt )
1970{
1971        int total;
1972        int nsite = 0;
1973        int atgcnum;
1974        int i, tmp;
1975        char *tmpseq, *tmpname;
1976        double atgcfreq;
1977        tmpname = AllocateCharVec( N );
1978        njob = countKUorWA( fp );
1979        searchKUorWA( fp );
1980        nlenmax = 0;
1981        *nlenminpt = 99999999;
1982        atgcnum = 0;
1983        total = 0;
1984        for( i=0; i<njob; i++ )
1985        {
1986                myfgets( tmpname, N-1, fp );
1987                tmpseq = load1SeqWithoutName_realloc( fp );
1988                tmp = countnogaplen( tmpseq );
1989                if( tmp > nlenmax ) nlenmax  = tmp;
1990                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
1991                atgcnum += countATGC( tmpseq, &nsite );
1992                total += nsite;
1993                free( tmpseq );
1994        }
1995        free( tmpname );
1996        atgcfreq = (double)atgcnum / total;
1997        fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1998        if( dorp == NOTSPECIFIED )
1999        {
2000                if( atgcfreq > 0.75 )   
2001                {
2002                        dorp = 'd';
2003                        upperCase = -1;
2004                }
2005                else                 
2006                {
2007                        dorp = 'p';
2008                        upperCase = 0;
2009                }
2010        }
2011}
2012
2013
2014void getnumlen_nogap_outallreg( FILE *fp, int *nlenminpt )
2015{
2016        int total;
2017        int nsite = 0;
2018        int atgcnum;
2019        int i, tmp;
2020        char *tmpseq, *tmpname;
2021        double atgcfreq;
2022        tmpname = AllocateCharVec( N );
2023        njob = countKUorWA( fp );
2024        searchKUorWA( fp );
2025        nlenmax = 0;
2026        *nlenminpt = 99999999;
2027        atgcnum = 0;
2028        total = 0;
2029        for( i=0; i<njob; i++ )
2030        {
2031                myfgets( tmpname, N-1, fp );
2032                fprintf( stdout, "%s\n", tmpname );
2033                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
2034                tmp = countnogaplen( tmpseq );
2035                fprintf( stdout, "%d\n", tmp );
2036                if( tmp > nlenmax ) nlenmax  = tmp;
2037                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2038                atgcnum += countATGC( tmpseq, &nsite );
2039                total += nsite;
2040                free( tmpseq );
2041        }
2042        free( tmpname );
2043        atgcfreq = (double)atgcnum / total;
2044        fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2045        if( dorp == NOTSPECIFIED )
2046        {
2047                if( atgcfreq > 0.75 )   
2048                {
2049                        dorp = 'd';
2050                        upperCase = -1;
2051                }
2052                else                 
2053                {
2054                        dorp = 'p';
2055                        upperCase = 0;
2056                }
2057        }
2058}
2059
2060static void escapehtml( char *res, char *ori, int maxlen )
2061{
2062        char *res0 = res;
2063        while( *ori )
2064        {
2065                if( *ori == '<' ) 
2066                {
2067                        strcpy( res, "&lt;" );
2068                        res += 3;
2069                }
2070                else if( *ori == '>' ) 
2071                {
2072                        strcpy( res, "&gt;" );
2073                        res += 3;
2074                }
2075                else if( *ori == '&' ) 
2076                {
2077                        strcpy( res, "&amp;" );
2078                        res += 4;
2079                }
2080                else if( *ori == '"' ) 
2081                {
2082                        strcpy( res, "&quot;" );
2083                        res += 5;
2084                }
2085                else if( *ori == ' ' ) 
2086                {
2087                        strcpy( res, "&nbsp;" );
2088                        res += 5;
2089                }
2090                else
2091                {
2092                        *res = *ori;
2093                }
2094                res++;
2095                ori++;
2096
2097                if( res - res0 -10 > N ) break;
2098        }
2099        *res = 0;
2100}
2101
2102void getnumlen_nogap_outallreg_web( FILE *fp, FILE *ofp, int *nlenminpt, int *isalignedpt )
2103{
2104        int total;
2105        int nsite = 0;
2106        int atgcnum;
2107        int alnlen = 0, alnlen_prev;
2108        int i, tmp, lennormalchar;
2109        char *tmpseq, *tmpname, *tmpname2;
2110        double atgcfreq;
2111        tmpname = AllocateCharVec( N );
2112        tmpname2 = AllocateCharVec( N );
2113        njob = countKUorWA( fp );
2114        searchKUorWA( fp );
2115        nlenmax = 0;
2116        *nlenminpt = 99999999;
2117        atgcnum = 0;
2118        total = 0;
2119        alnlen_prev = -1;
2120        *isalignedpt = 1;
2121        for( i=0; i<njob; i++ )
2122        {
2123                myfgets( tmpname, N-1, fp );
2124                tmpname2[0] = tmpname[0];
2125                escapehtml( tmpname2+1, tmpname+1, N );
2126//              fprintf( stdout, "%s\n", tmpname );
2127//              fprintf( stdout, "%s\n", tmpname2 );
2128//              exit(1);
2129                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
2130                tmp = countnogaplen( tmpseq );
2131//              fprintf( stdout, "%d\n", tmp );
2132                if( tmp > nlenmax ) nlenmax  = tmp;
2133                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2134                atgcnum += countATGC( tmpseq, &nsite );
2135                total += nsite;
2136
2137                alnlen = strlen( tmpseq );
2138//              fprintf( stdout, "##### alnlen, alnlen_prev = %d, %d\n", alnlen, alnlen_prev );
2139                if( i>0 && alnlen_prev != alnlen ) *isalignedpt = 0;
2140                alnlen_prev = alnlen;
2141
2142                atgcfreq = (double)atgcnum / total;
2143//              fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2144//              if( dorp == NOTSPECIFIED ) // you kentou
2145                {
2146                        if( atgcfreq > 0.75 )   
2147                        {
2148                                dorp = 'd';
2149                                upperCase = -1;
2150                        }
2151                        else                 
2152                        {
2153                                dorp = 'p';
2154                                upperCase = 0;
2155                        }
2156                }
2157       
2158                if( dorp == 'd' ) lennormalchar = countnormalletters( tmpseq, "atgcuATGCU" );
2159                else              lennormalchar = countnormalletters( tmpseq, "ARNDCQEGHILKMFPSTWYVarndcqeghilkmfpstwyv" );
2160                free( tmpseq );
2161
2162                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 );
2163                fprintf( ofp, "<span id='t%d-0' style='display:none'>", i );
2164                fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"t%d\")'>+reg</a>", i );
2165                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 );
2166                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 );
2167//              fprintf( ofp, "  Sequence Length:<input type='text' name='l%d' size='8' value='%d' readonly='readonly'>", i, tmp );
2168                fprintf( ofp, "\n</span>" );
2169                fprintf( ofp, "<span id='t%d-1' style='display:none'>", i );
2170                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 );
2171                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 );
2172                fprintf( ofp, "\n</span>" );
2173                fprintf( ofp, "<span id='t%d-2' style='display:none'>", i );
2174                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 );
2175                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 );
2176                fprintf( ofp, "\n</span>" );
2177                fprintf( ofp, "<span id='t%d-3' style='display:none'>", i );
2178                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 );
2179                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 );
2180                fprintf( ofp, "\n</span>" );
2181                fprintf( ofp, "<span id='t%d-4' style='display:none'>", i );
2182                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 );
2183                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 );
2184                fprintf( ofp, "\n</span>" );
2185        }
2186        free( tmpname );
2187        free( tmpname2 );
2188        atgcfreq = (double)atgcnum / total;
2189        fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2190//      if( dorp == NOTSPECIFIED ) // you kentou
2191        {
2192                if( atgcfreq > 0.75 )   
2193                {
2194                        dorp = 'd';
2195                        upperCase = -1;
2196                }
2197                else                 
2198                {
2199                        dorp = 'p';
2200                        upperCase = 0;
2201                }
2202        }
2203        fprintf( ofp, "\n" );
2204        if( *isalignedpt )
2205        {
2206                fprintf( ofp, "<span id='tall-0' style='display:none'>" );
2207                fprintf( ofp, "Cut the alignment\n" );
2208                fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"tall\")'>+reg</a>" );
2209                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 );
2210                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-0'><input type='checkbox' name='rall-0' id='rall-0'>Reverse</label>" );
2211                fprintf( ofp, "  Alignment length:<input type='text' name='lall' size='8' value='%d' readonly='readonly'>", alnlen );
2212                fprintf( ofp, "\n</span>" );
2213                fprintf( ofp, "<span id='tall-1' style='display:none'>" );
2214                fprintf( ofp, "      Begin:<input type='text' name='ball-1' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2215                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-1'><input type='checkbox' name='rall-1' id='rall-1'>Reverse</label>" );
2216                fprintf( ofp, "\n</span>" );
2217                fprintf( ofp, "<span id='tall-2' style='display:none'>" );
2218                fprintf( ofp, "      Begin:<input type='text' name='ball-2' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2219                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-2'><input type='checkbox' name='rall-2' id='rall-2'>Reverse</label>" );
2220                fprintf( ofp, "\n</span>" );
2221                fprintf( ofp, "<span id='tall-3' style='display:none'>" );
2222                fprintf( ofp, "      Begin:<input type='text' name='ball-3' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2223                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-3'><input type='checkbox' name='rall-3' id='rall-3'>Reverse</label>" );
2224                fprintf( ofp, "\n</span>" );
2225                fprintf( ofp, "<span id='tall-4' style='display:none'>" );
2226                fprintf( ofp, "      Begin:<input type='text' name='ball-4' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2227                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-4'><input type='checkbox' name='rall-4' id='rall-4'>Reverse</label>" );
2228                fprintf( ofp, "\n</span>" );
2229        }
2230
2231}
2232
2233void getnumlen( FILE *fp )
2234{
2235        int total;
2236        int nsite = 0;
2237        int atgcnum;
2238        int i, tmp;
2239        char *tmpseq;
2240        char *tmpname;
2241        double atgcfreq;
2242        tmpname = AllocateCharVec( N );
2243        njob = countKUorWA( fp );
2244        searchKUorWA( fp );
2245        nlenmax = 0;
2246        atgcnum = 0;
2247        total = 0;
2248        for( i=0; i<njob; i++ )
2249        {
2250                myfgets( tmpname, N-1, fp );
2251                tmpseq = load1SeqWithoutName_realloc( fp );
2252                tmp = strlen( tmpseq );
2253                if( tmp > nlenmax ) nlenmax  = tmp;
2254                atgcnum += countATGC( tmpseq, &nsite );
2255                total += nsite;
2256//              fprintf( stderr, "##### total = %d\n", total );
2257                free( tmpseq );
2258        }
2259
2260
2261        atgcfreq = (double)atgcnum / total;
2262//      fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2263        if( dorp == NOTSPECIFIED )
2264        {
2265                if( atgcfreq > 0.75 )   
2266                {
2267                        dorp = 'd';
2268                        upperCase = -1;
2269                }
2270                else                 
2271                {
2272                        dorp = 'p';
2273                        upperCase = 0;
2274                }
2275        }
2276        free( tmpname );
2277}
2278       
2279
2280
2281void WriteGapFill( FILE *fp, int locnjob, char name[][B], int nlen[M], char **aseq )
2282{
2283        static char b[N];
2284        int i, j;
2285        int nalen[M];
2286        static char gap[N];
2287        static char buff[N];
2288
2289#if IODEBUG
2290        fprintf( stderr, "IMAKARA KAKU\n" );
2291#endif
2292        nlenmax = 0;
2293        for( i=0; i<locnjob; i++ )
2294        {
2295                int len = strlen( aseq[i] );
2296                if( nlenmax < len ) nlenmax = len;
2297        }
2298
2299        for( i=0; i<nlenmax; i++ ) gap[i] = '-';
2300        gap[nlenmax] = 0;
2301
2302        fprintf( fp, "%5d", locnjob );
2303        fprintf( fp, "\n" );
2304
2305        for( i=0; i<locnjob; i++ )
2306        {
2307                strcpy( buff, aseq[i] );
2308                strncat( buff, gap, nlenmax-strlen( aseq[i] ) );
2309                buff[nlenmax] = 0;
2310                nalen[i] = strlen( buff );
2311                fprintf( fp, "%s\n", name[i] );
2312                fprintf( fp, "%5d\n", nalen[i] );
2313                for( j=0; j<nalen[i]; j=j+C )
2314                {
2315                        strncpy_caseC( b, buff+j, C ); b[C] = 0;
2316                        fprintf( fp, "%s\n",b );
2317                }
2318        }
2319#if DEBUG
2320        fprintf( stderr, "nalen[0] = %d\n", nalen[0] );
2321#endif
2322#if IODEBUG
2323        fprintf( stderr, "KAKIOWATTA\n" );
2324#endif
2325}
2326
2327void writeDataforgaln( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2328{
2329        int i, j;
2330        int nalen;
2331
2332        for( i=0; i<locnjob; i++ )
2333        {
2334                nalen = strlen( aseq[i] );
2335                fprintf( fp, ">%s\n", name[i]+1 );
2336                for( j=0; j<nalen; j=j+C )
2337                {
2338#if 0
2339                        strncpy( b, aseq[i]+j, C ); b[C] = 0;
2340                        fprintf( fp, "%s\n",b );
2341#else
2342                        fprintf( fp, "%.*s\n", C, aseq[i]+j );
2343#endif
2344                }
2345        }
2346}
2347
2348void writeData_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2349{
2350        int i, j;
2351        int nalen;
2352
2353        for( i=0; i<locnjob; i++ )
2354        {
2355#if DEBUG
2356                fprintf( stderr, "i = %d in writeData\n", i );
2357#endif
2358                nalen = strlen( aseq[i] );
2359                fprintf( fp, ">%s\n", name[i]+1 );
2360                for( j=0; j<nalen; j=j+C )
2361                {
2362#if 0
2363                        strncpy( b, aseq[i]+j, C ); b[C] = 0;
2364                        fprintf( fp, "%s\n",b );
2365#else
2366                        fprintf( fp, "%.*s\n", C, aseq[i]+j );
2367#endif
2368                }
2369        }
2370}
2371
2372void writeData( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq )
2373{
2374        int i, j;
2375        int nalen;
2376
2377        for( i=0; i<locnjob; i++ )
2378        {
2379#if DEBUG
2380                fprintf( stderr, "i = %d in writeData\n", i );
2381#endif
2382                nalen = strlen( aseq[i] );
2383                fprintf( fp, ">%s\n", name[i]+1 );
2384                for( j=0; j<nalen; j=j+C )
2385                {
2386#if 0
2387                        strncpy( b, aseq[i]+j, C ); b[C] = 0;
2388                        fprintf( fp, "%s\n",b );
2389#else
2390                        fprintf( fp, "%.*s\n", C, aseq[i]+j );
2391#endif
2392                }
2393        }
2394}
2395
2396
2397void write1seq( FILE *fp, char *aseq )
2398{
2399        int j;
2400        int nalen;
2401
2402        nalen = strlen( aseq );
2403        for( j=0; j<nalen; j=j+C )
2404                fprintf( fp, "%.*s\n", C, aseq+j );
2405}
2406
2407void readhat2_floathalf_part_pointer( FILE *fp, int nseq, int nadd, char **name, float **mtx )
2408{
2409    int i, j, nseq0, norg;
2410    char b[B];
2411
2412    fgets( b, B, fp );
2413    fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) 
2414        {
2415                fprintf( stderr, "%d != %d\n", nseq, nseq0 );
2416                ErrorExit( "hat2 is wrong." );
2417        }
2418    fgets( b, B, fp );
2419    for( i=0; i<nseq; i++ )
2420    {
2421#if 0
2422        getaline_fp_eof( b, B, fp );
2423#else
2424                myfgets( b, B-2, fp );
2425#endif
2426#if 0
2427                j = MIN( strlen( b+6 ), 10 );
2428        if( strncmp( name[i], b+6 , j ) )
2429                {
2430                        fprintf( stderr, "Error in hat2\n" );
2431                        fprintf( stderr, "%s != %s\n", b, name[i] );
2432                        exit( 1 );
2433                }
2434#endif
2435    }
2436        norg = nseq-nadd;
2437    for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ )
2438    {
2439        mtx[i][j] = ( input_new( fp, D ) );
2440    }
2441}
2442
2443void readhat2_floathalf_pointer( FILE *fp, int nseq, char **name, float **mtx )
2444{
2445    int i, j, nseq0;
2446    char b[B];
2447
2448    fgets( b, B, fp );
2449    fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) 
2450        {
2451                fprintf( stderr, "%d != %d\n", nseq, nseq0 );
2452                ErrorExit( "hat2 is wrong." );
2453        }
2454    fgets( b, B, fp );
2455    for( i=0; i<nseq; i++ )
2456    {
2457#if 0
2458        getaline_fp_eof( b, B, fp );
2459#else
2460                myfgets( b, B-2, fp );
2461#endif
2462#if 0
2463                j = MIN( strlen( b+6 ), 10 );
2464        if( strncmp( name[i], b+6 , j ) )
2465                {
2466                        fprintf( stderr, "Error in hat2\n" );
2467                        fprintf( stderr, "%s != %s\n", b, name[i] );
2468                        exit( 1 );
2469                }
2470#endif
2471    }
2472    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2473    {
2474        mtx[i][j-i] = ( input_new( fp, D ) );
2475    }
2476}
2477void readhat2_floathalf( FILE *fp, int nseq, char name[M][B], float **mtx )
2478{
2479    int i, j, nseq0;
2480    char b[B];
2481
2482    fgets( b, B, fp );
2483    fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2484    fgets( b, B, fp );
2485    for( i=0; i<nseq; i++ )
2486    {
2487#if 0
2488        getaline_fp_eof( b, B, fp );
2489#else
2490                myfgets( b, B-2, fp );
2491#endif
2492#if 0
2493                j = MIN( strlen( b+6 ), 10 );
2494        if( strncmp( name[i], b+6 , j ) )
2495                {
2496                        fprintf( stderr, "Error in hat2\n" );
2497                        fprintf( stderr, "%s != %s\n", b, name[i] );
2498                        exit( 1 );
2499                }
2500#endif
2501    }
2502    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2503    {
2504        mtx[i][j-i] = ( input_new( fp, D ) );
2505    }
2506}
2507void readhat2_float( FILE *fp, int nseq, char name[M][B], float **mtx )
2508{
2509    int i, j, nseq0;
2510    char b[B];
2511
2512    fgets( b, B, fp );
2513    fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2514    fgets( b, B, fp );
2515    for( i=0; i<nseq; i++ )
2516    {
2517#if 0
2518        getaline_fp_eof( b, B, fp );
2519#else
2520                myfgets( b, B-2, fp );
2521#endif
2522#if 0
2523                j = MIN( strlen( b+6 ), 10 );
2524        if( strncmp( name[i], b+6 , j ) )
2525                {
2526                        fprintf( stderr, "Error in hat2\n" );
2527                        fprintf( stderr, "%s != %s\n", b, name[i] );
2528                        exit( 1 );
2529                }
2530#endif
2531    }
2532    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2533    {
2534        mtx[i][j] = ( input_new( fp, D ) );
2535    }
2536}
2537void readhat2_int( FILE *fp, int nseq, char name[M][B], int **mtx )
2538{
2539    int i, j, nseq0;
2540    char b[B];
2541
2542    fgets( b, B, fp );
2543    fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2544    fgets( b, B, fp );
2545    for( i=0; i<nseq; i++ )
2546    {
2547#if 0
2548        getaline_fp_eof( b, B, fp );
2549#else
2550                myfgets( b, B-2, fp );
2551#endif
2552#if 0
2553                j = MIN( strlen( b+6 ), 10 );
2554        if( strncmp( name[i], b+6 , j ) )
2555                {
2556                        fprintf( stderr, "Error in hat2\n" );
2557                        fprintf( stderr, "%s != %s\n", b, name[i] );
2558                        exit( 1 );
2559                }
2560#endif
2561    }
2562    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2563    {
2564        mtx[i][j] = (int)( input_new( fp, D ) * INTMTXSCALE + 0.5 );
2565    }
2566}
2567
2568void readhat2_pointer( FILE *fp, int nseq, char **name, double **mtx )
2569{
2570    int i, j, nseq0;
2571    char b[B];
2572
2573    fgets( b, B, fp );
2574    fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2575    fgets( b, B, fp );
2576    for( i=0; i<nseq; i++ )
2577    {
2578#if 0
2579        getaline_fp_eof( b, B, fp );
2580#else
2581                myfgets( b, B-2, fp );
2582#endif
2583#if 0
2584                j = MIN( strlen( b+6 ), 10 );
2585        if( strncmp( name[i], b+6 , j ) )
2586                {
2587                        fprintf( stderr, "Error in hat2\n" );
2588                        fprintf( stderr, "%s != %s\n", b, name[i] );
2589                        exit( 1 );
2590                }
2591#endif
2592    }
2593    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2594    {
2595        mtx[i][j] = (double)input_new( fp, D);
2596    }
2597}
2598void readhat2( FILE *fp, int nseq, char name[M][B], double **mtx )
2599{
2600    int i, j, nseq0;
2601    char b[B];
2602
2603    fgets( b, B, fp );
2604    fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2605    fgets( b, B, fp );
2606    for( i=0; i<nseq; i++ )
2607    {
2608#if 0
2609        getaline_fp_eof( b, B, fp );
2610#else
2611                myfgets( b, B-2, fp );
2612#endif
2613#if 0
2614                j = MIN( strlen( b+6 ), 10 );
2615        if( strncmp( name[i], b+6 , j ) )
2616                {
2617                        fprintf( stderr, "Error in hat2\n" );
2618                        fprintf( stderr, "%s != %s\n", b, name[i] );
2619                        exit( 1 );
2620                }
2621#endif
2622    }
2623    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2624    {
2625        mtx[i][j] = (double)input_new( fp, D);
2626    }
2627}
2628
2629void WriteFloatHat2_pointer_halfmtx( FILE *hat2p, int locnjob, char **name, float **mtx )
2630{
2631        int i, j, ijsa;
2632        double max = 0.0;
2633        for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2634
2635        fprintf( hat2p, "%5d\n", 1 );
2636        fprintf( hat2p, "%5d\n", locnjob );
2637        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2638
2639        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2640        for( i=0; i<locnjob; i++ )
2641        {
2642                for( j=i+1; j<njob; j++ )
2643                {
2644                        fprintf( hat2p, "%#6.3f", mtx[i][j-i] );
2645                        ijsa = j-i;
2646                        if( ijsa % 12 == 0 || ijsa == locnjob-i-1 ) fprintf( hat2p, "\n" );
2647                }
2648        }
2649}
2650
2651void WriteFloatHat2_pointer( FILE *hat2p, int locnjob, char **name, float **mtx )
2652{
2653        int i, j;
2654        double max = 0.0;
2655        for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2656
2657        fprintf( hat2p, "%5d\n", 1 );
2658        fprintf( hat2p, "%5d\n", locnjob );
2659        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2660
2661        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2662        for( i=0; i<locnjob; i++ )
2663        {
2664                for( j=1; j<locnjob-i; j++ ) 
2665                {
2666                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2667                        if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2668                }
2669        }
2670}
2671
2672void WriteFloatHat2( FILE *hat2p, int locnjob, char name[M][B], float **mtx )
2673{
2674        int i, j;
2675        double max = 0.0;
2676        for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2677
2678        fprintf( hat2p, "%5d\n", 1 );
2679        fprintf( hat2p, "%5d\n", locnjob );
2680        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2681
2682        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2683        for( i=0; i<locnjob; i++ )
2684        {
2685                for( j=1; j<locnjob-i; j++ ) 
2686                {
2687                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2688                        if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2689                }
2690        }
2691}
2692
2693void WriteHat2_int( FILE *hat2p, int locnjob, char name[M][B], int **mtx )
2694{
2695        int i, j;
2696        double max = 0.0;
2697        for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2698        max /= INTMTXSCALE;
2699
2700        fprintf( hat2p, "%5d\n", 1 );
2701        fprintf( hat2p, "%5d\n", locnjob );
2702        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2703
2704        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2705        for( i=0; i<locnjob-1; i++ )
2706        {
2707                for( j=i+1; j<locnjob; j++ ) 
2708                {
2709                        fprintf( hat2p, "%#6.3f", (float)mtx[i][j] / INTMTXSCALE );
2710                        if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2711                }
2712        }
2713}
2714
2715void WriteHat2_part_pointer( FILE *hat2p, int locnjob, int nadd, char **name, double **mtx )
2716{
2717        int i, j;
2718        int norg = locnjob-nadd;
2719        double max = 0.0;
2720//      for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2721
2722        fprintf( hat2p, "%5d\n", 1 );
2723        fprintf( hat2p, "%5d\n", locnjob );
2724        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2725
2726        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2727        for( i=0; i<norg; i++ )
2728        {
2729                for( j=0; j<nadd; j++ ) 
2730                {
2731                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2732                        if( (j+1) % 12 == 0 || j == nadd-1 ) fprintf( hat2p, "\n" );
2733                }
2734        }
2735}
2736
2737void WriteHat2_pointer( FILE *hat2p, int locnjob, char **name, double **mtx )
2738{
2739        int i, j;
2740        double max = 0.0;
2741        for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2742
2743        fprintf( hat2p, "%5d\n", 1 );
2744        fprintf( hat2p, "%5d\n", locnjob );
2745        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2746
2747        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2748        for( i=0; i<locnjob-1; i++ )
2749        {
2750                for( j=i+1; j<locnjob; j++ ) 
2751                {
2752                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2753                        if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2754                }
2755        }
2756}
2757
2758void WriteHat2( FILE *hat2p, int locnjob, char name[M][B], double **mtx )
2759{
2760        int i, j;
2761        double max = 0.0;
2762        for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2763
2764        fprintf( hat2p, "%5d\n", 1 );
2765        fprintf( hat2p, "%5d\n", locnjob );
2766        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2767
2768        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2769        for( i=0; i<locnjob-1; i++ )
2770        {
2771                for( j=i+1; j<locnjob; j++ ) 
2772                {
2773                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2774                        if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2775                }
2776        }
2777}
2778
2779#if 0
2780void WriteHat2plain( FILE *hat2p, int locnjob, double **mtx )
2781{
2782        int i, j, ilim;
2783
2784        ilim = locnjob-1;
2785        for( i=0; i<ilim; i++ )
2786        {
2787                fprintf( hat2p, "%d-%d d=%.3f\n", i+1, i+1, 0.0 );
2788                for( j=i+1; j<locnjob; j++ )
2789                {
2790                        fprintf( hat2p, "%d-%d d=%.3f\n", i+1, j+1, mtx[i][j] );
2791                }
2792        }
2793}
2794#endif
2795
2796int ReadFasta_sub( FILE *fp, double *dis, int nseq, char name[M][B] )
2797{
2798    int i, count=0;
2799    char b[B];
2800    int junban[M];
2801
2802    count = 0;
2803    for( i=0; i<10000000 && count<nseq; i++ )
2804    {
2805        fgets( b, B-1, fp );
2806        if( !strncmp( "+==========+", b, 12 ) )
2807        {
2808            junban[count] = atoi( b+12 );
2809            count++;
2810        }
2811    }
2812
2813        for( i=0; i<nseq; i++ ) dis[i] = 0.0;
2814    count = 0;
2815    for( i=0; i<100000 && count<nseq; i++ )
2816    {
2817                if( fgets( b, B-1, fp ) ) break;
2818        if( !strncmp( name[junban[count]], b, 20  ) )
2819        {
2820            fgets( b, B-1, fp );
2821            dis[junban[count]] = atof( b );
2822            count++;
2823        }
2824    }
2825    return 0;
2826}
2827
2828
2829int ReadSsearch( FILE *fp, double *dis, int nseq, char name[M][B] )
2830{
2831    int i, count=0;
2832    char b[B];
2833    int junban[M];
2834        int opt;
2835
2836    count = 0;
2837    for( i=0; i<10000000 && count<nseq; i++ )
2838    {
2839        fgets( b, B-1, fp );
2840        if( !strncmp( "+==========+", b, 12 ) )
2841        {
2842            junban[count] = atoi( b+12 );
2843                        sscanf( b+75, "%d", &opt ); 
2844            dis[junban[count]] = (double)opt;
2845            count++;
2846        }
2847    }
2848
2849/*
2850    count = 0;
2851    for( i=0; i<100000 && count<nseq; i++ )
2852    {
2853        fgets( b, B-1, fp );
2854        if( !strncmp( name[junban[count]], b, 20  ) )
2855        {
2856            dis[junban[count]] = atof( b+65 );
2857            count++;
2858        }
2859    }
2860*/
2861    return 0;
2862}
2863
2864int ReadBlastm7_avscore( FILE *fp, double *dis, int nin )
2865{
2866    int count=0;
2867    char b[B];
2868        char *pt;
2869    int *junban;
2870        double score, sumscore;
2871        double len, sumlen;
2872        int qstart, qend, tstart, tend;
2873        double scorepersite;
2874        static char qal[N], tal[N], al[N];
2875        int nlocalhom;
2876
2877        junban = calloc( nin, sizeof( int ) );
2878
2879        count = 0;
2880        sumscore = 0.0;
2881        sumlen = 0.0;
2882        score = 0.0;
2883        len = 0.0;
2884        scorepersite = 0.0; // by D.Mathog, a guess
2885    while( 1 )
2886        {
2887
2888                if( feof( fp ) ) break;
2889
2890                while( fgets( b, B-1, fp ) )
2891                {
2892                        if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
2893                }
2894
2895                if( !strncmp( "          <Hit_def>", b, 19 ) )
2896                {
2897                        junban[count] = atoi( b+31 );
2898                        nlocalhom = 0;
2899                }
2900
2901
2902                while( fgets( b, B-1, fp ) )
2903                        if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
2904                pt = b + 25;
2905                score = atof( pt );
2906                sumscore += score;
2907
2908
2909                while( fgets( b, B-1, fp ) )
2910                        if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
2911                pt = b + 30;
2912                qstart = atoi( pt ) - 1;
2913
2914
2915                while( fgets( b, B-1, fp ) )
2916                        if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
2917                pt = b + 28;
2918                qend = atoi( pt ) - 1;
2919
2920
2921                while( fgets( b, B-1, fp ) )
2922                        if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
2923                pt = b + 28;
2924                tstart = atoi( pt ) - 1;
2925
2926
2927                while( fgets( b, B-1, fp ) )
2928                        if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
2929                pt = b + 26;
2930                tend = atoi( pt ) - 1;
2931
2932
2933                while( fgets( b, B-1, fp ) )
2934                        if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
2935                pt = b + 29;
2936                len = atoi( pt );
2937                sumlen += len;
2938
2939
2940                while( fgets( al, N-100, fp ) )
2941                        if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
2942
2943                strcpy( qal, al+24 );
2944                pt = qal;
2945                while( *++pt != '<' )
2946                        ;
2947                *pt = 0;
2948
2949
2950                while( fgets( al, N-100, fp ) )
2951                        if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
2952
2953                strcpy( tal, al+24 );
2954                pt = tal;
2955                while( *++pt != '<' )
2956                        ;
2957                *pt = 0;
2958
2959
2960//              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 );
2961
2962
2963                while( fgets( b, B-1, fp ) )
2964                        if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
2965
2966
2967                fgets( b, B-1, fp );
2968
2969
2970                if( !strncmp( "          </Hit_hsps>", b, 21 ) )
2971                {
2972                        dis[junban[count++]] = sumscore;
2973                        sumscore = 0.0;
2974                        fgets( b, B-1, fp );
2975                        fgets( b, B-1, fp );
2976                        scorepersite = sumscore / sumlen;
2977                        if( scorepersite != (int)scorepersite )
2978                        {
2979                                fprintf( stderr, "ERROR! sumscore=%f, sumlen=%f, and scorepersite=%f\n", sumscore, sumlen, scorepersite );
2980                                exit( 1 );
2981                        }
2982
2983                        if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
2984                }
2985        }
2986
2987        free( junban );
2988
2989    return (int)scorepersite;
2990}
2991int ReadBlastm7_scoreonly( FILE *fp, double *dis, int nin )
2992{
2993    int count=0;
2994    char b[B];
2995        char *pt;
2996    int *junban;
2997        int overlapaa;
2998        double score, sumscore;
2999        int qstart, qend, tstart, tend;
3000        static char qal[N], tal[N], al[N];
3001        int nlocalhom;
3002
3003        junban = calloc( nin, sizeof( int ) );
3004
3005        count = 0;
3006        sumscore = 0.0;
3007        score = 0.0;
3008    while( 1 )
3009        {
3010
3011                if( feof( fp ) ) break;
3012
3013                while( fgets( b, B-1, fp ) )
3014                {
3015                        if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
3016                }
3017
3018                if( !strncmp( "          <Hit_def>", b, 19 ) )
3019                {
3020                        junban[count] = atoi( b+31 );
3021                        nlocalhom = 0;
3022                }
3023
3024
3025                while( fgets( b, B-1, fp ) )
3026                        if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
3027                pt = b + 25;
3028                score = atof( pt );
3029                sumscore += score;
3030
3031
3032                while( fgets( b, B-1, fp ) )
3033                        if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
3034                pt = b + 30;
3035                qstart = atoi( pt ) - 1;
3036
3037
3038                while( fgets( b, B-1, fp ) )
3039                        if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
3040                pt = b + 28;
3041                qend = atoi( pt ) - 1;
3042
3043
3044                while( fgets( b, B-1, fp ) )
3045                        if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
3046                pt = b + 28;
3047                tstart = atoi( pt ) - 1;
3048
3049
3050                while( fgets( b, B-1, fp ) )
3051                        if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
3052                pt = b + 26;
3053                tend = atoi( pt ) - 1;
3054
3055
3056                while( fgets( b, B-1, fp ) )
3057                        if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
3058                pt = b + 29;
3059                overlapaa = atoi( pt );
3060
3061
3062                while( fgets( al, N-100, fp ) )
3063                        if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
3064
3065                strcpy( qal, al+24 );
3066                pt = qal;
3067                while( *++pt != '<' )
3068                        ;
3069                *pt = 0;
3070
3071
3072                while( fgets( al, N-100, fp ) )
3073                        if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
3074
3075                strcpy( tal, al+24 );
3076                pt = tal;
3077                while( *++pt != '<' )
3078                        ;
3079                *pt = 0;
3080
3081
3082//              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 );
3083
3084//              nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
3085
3086                while( fgets( b, B-1, fp ) )
3087                        if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
3088
3089
3090                fgets( b, B-1, fp );
3091
3092
3093                if( !strncmp( "          </Hit_hsps>", b, 21 ) )
3094                {
3095                        dis[junban[count++]] = sumscore;
3096                        sumscore = 0.0;
3097                        fgets( b, B-1, fp );
3098                        fgets( b, B-1, fp );
3099                        if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
3100                }
3101        }
3102
3103        free( junban );
3104
3105    return count;
3106}
3107
3108int ReadBlastm7( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3109{
3110    int count=0;
3111    char b[B];
3112        char *pt;
3113    static int junban[M];
3114        int overlapaa;
3115        double score, sumscore;
3116        int qstart, qend, tstart, tend;
3117        static char qal[N], tal[N], al[N];
3118        int nlocalhom;
3119
3120
3121
3122        count = 0;
3123        sumscore = 0.0;
3124        score = 0.0;
3125        nlocalhom = 0;
3126    while( 1 )
3127        {
3128
3129                if( feof( fp ) ) break;
3130
3131                while( fgets( b, B-1, fp ) )
3132                {
3133                        if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
3134                }
3135
3136                if( !strncmp( "          <Hit_def>", b, 19 ) )
3137                {
3138                        junban[count] = atoi( b+31 );
3139                        nlocalhom = 0;
3140                }
3141
3142
3143                while( fgets( b, B-1, fp ) )
3144                        if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
3145                pt = b + 25;
3146                score = atof( pt );
3147                sumscore += score;
3148
3149
3150                while( fgets( b, B-1, fp ) )
3151                        if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
3152                pt = b + 30;
3153                qstart = atoi( pt ) - 1;
3154
3155
3156                while( fgets( b, B-1, fp ) )
3157                        if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
3158                pt = b + 28;
3159                qend = atoi( pt ) - 1;
3160
3161
3162                while( fgets( b, B-1, fp ) )
3163                        if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
3164                pt = b + 28;
3165                tstart = atoi( pt ) - 1;
3166
3167
3168                while( fgets( b, B-1, fp ) )
3169                        if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
3170                pt = b + 26;
3171                tend = atoi( pt ) - 1;
3172
3173
3174                while( fgets( b, B-1, fp ) )
3175                        if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
3176                pt = b + 29;
3177                overlapaa = atoi( pt );
3178
3179
3180                while( fgets( al, N-100, fp ) )
3181                        if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
3182
3183                strcpy( qal, al+24 );
3184                pt = qal;
3185                while( *++pt != '<' )
3186                        ;
3187                *pt = 0;
3188
3189
3190                while( fgets( al, N-100, fp ) )
3191                        if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
3192
3193                strcpy( tal, al+24 );
3194                pt = tal;
3195                while( *++pt != '<' )
3196                        ;
3197                *pt = 0;
3198
3199
3200//              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 );
3201
3202                nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
3203
3204                while( fgets( b, B-1, fp ) )
3205                        if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
3206
3207
3208                fgets( b, B-1, fp );
3209
3210
3211                if( !strncmp( "          </Hit_hsps>", b, 21 ) )
3212                {
3213                        dis[junban[count++]] = sumscore;
3214                        sumscore = 0.0;
3215                        fgets( b, B-1, fp );
3216                        fgets( b, B-1, fp );
3217                        if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
3218                }
3219        }
3220    return count;
3221}
3222
3223int ReadFasta34noalign( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3224{
3225    int count=0;
3226    char b[B];
3227        char *pt;
3228    static int junban[M];
3229        int opt;
3230        double z, bits;
3231
3232
3233    count = 0;
3234#if 0
3235    for( i=0; i<10000000 && count<nseq; i++ )
3236#else
3237    while( !feof( fp ) )
3238#endif
3239    {
3240        fgets( b, B-1, fp );
3241        if( !strncmp( "+==========+", b, 12 ) )
3242        {
3243            junban[count] = atoi( b+12 );
3244
3245                        pt = strchr( b, ')' ) + 1;
3246                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3247            dis[junban[count]] = (double)opt;
3248            count++;
3249
3250        }
3251    }
3252
3253    return count;
3254}
3255int ReadFasta34m10_nuc( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3256{
3257    int count=0;
3258    char b[B];
3259        char *pt;
3260    static int junban[M];
3261        int overlapaa;
3262        int opt, qstart, qend, tstart, tend;
3263        double z, bits;
3264        int qal_display_start, tal_display_start;
3265        static char qal[N], tal[N];
3266        char *qal2, *tal2;
3267        int c;
3268
3269
3270    count = 0;
3271#if 0
3272    for( i=0; i<10000000 && count<nseq; i++ )
3273#else
3274    while( !feof( fp ) )
3275#endif
3276    {
3277        fgets( b, B-1, fp );
3278        if( !strncmp( "+==========+", b, 12 ) )
3279        {
3280            junban[count] = atoi( b+12 );
3281
3282                        if( strchr( b, 'r' ) ) continue;
3283
3284                        pt = strchr( b, ']' ) + 1;
3285                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3286            dis[junban[count]] = (double)opt;
3287            count++;
3288
3289        }
3290                else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3291                {
3292                        break;
3293                }
3294
3295    }
3296        if( !count ) return -1;
3297
3298        count = 0;
3299    while( 1 )
3300        {
3301                if( strncmp( ">>+==========+", b, 14 ) )
3302                {
3303                        fgets( b, B-1, fp );
3304                        if( feof( fp ) ) break;
3305                        continue;
3306                }
3307                junban[count++] = atoi( b+14 );
3308//              fprintf( stderr, "t = %d\n", atoi( b+14 ) );
3309                while( fgets( b, B-1, fp ) )
3310                        if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
3311                pt = strstr( b, ":" ) +1;
3312                opt = atoi( pt );
3313
3314
3315                while( fgets( b, B-1, fp ) )
3316                        if( !strncmp( "_overlap:", b+4, 9 ) ) break;
3317                pt = strstr( b, ":" ) +1;
3318                overlapaa = atoi( pt );
3319
3320                while( fgets( b, B-1, fp ) )
3321                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3322                pt = strstr( b, ":" ) +1;
3323                qstart = atoi( pt ) - 1;
3324
3325                while( fgets( b, B-1, fp ) )
3326                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3327                pt = strstr( b, ":" ) +1;
3328                qend = atoi( pt ) - 1;
3329
3330                while( fgets( b, B-1, fp ) )
3331                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3332                pt = strstr( b, ":" ) +1;
3333                qal_display_start = atoi( pt ) - 1;
3334
3335                pt = qal;
3336                while( (c = fgetc( fp )) )
3337                {
3338                        if( c == '>' ) 
3339                        {
3340                                ungetc( c, fp );
3341                                break;
3342                        }
3343                        if( isalpha( c ) || c == '-' ) 
3344                        *pt++ = c;
3345                }
3346                *pt = 0;
3347
3348                while( fgets( b, B-1, fp ) )
3349                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3350                pt = strstr( b, ":" ) + 1;
3351                tstart = atoi( pt ) - 1;
3352
3353                while( fgets( b, B-1, fp ) )
3354                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3355                pt = strstr( b, ":" ) + 1;
3356                tend = atoi( pt ) - 1;
3357
3358                while( fgets( b, B-1, fp ) )
3359                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3360                pt = strstr( b, ":" ) + 1;
3361                tal_display_start = atoi( pt ) - 1;
3362
3363                pt = tal;
3364                while( ( c = fgetc( fp ) ) )
3365                {
3366                        if( c == '>' ) 
3367                        {
3368                                ungetc( c, fp );
3369                                break;
3370                        }
3371                        if( isalpha( c ) || c == '-' ) 
3372                        *pt++ = c;
3373                }
3374                *pt = 0;
3375
3376//              fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
3377//              fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
3378
3379//              fprintf( stderr, "qal = %s\n", qal );
3380//              fprintf( stderr, "tal = %s\n", tal );
3381
3382                qal2 = cutal( qal, qal_display_start, qstart, qend );
3383                tal2 = cutal( tal, tal_display_start, tstart, tend );
3384
3385//              fprintf( stderr, "qal2 = %s\n", qal2 );
3386//              fprintf( stderr, "tal2 = %s\n", tal2 );
3387
3388//              fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
3389                putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
3390        }
3391//      fprintf( stderr, "count = %d\n", count );
3392    return count;
3393}
3394int ReadFasta34m10( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3395{
3396    int count=0;
3397    char b[B];
3398        char *pt;
3399    static int junban[M];
3400        int overlapaa;
3401        int opt, qstart, qend, tstart, tend;
3402        double z, bits;
3403        int qal_display_start, tal_display_start;
3404        static char qal[N], tal[N];
3405        char *qal2, *tal2;
3406        int c;
3407
3408
3409    count = 0;
3410#if 0
3411    for( i=0; i<10000000 && count<nseq; i++ )
3412#else
3413    while( !feof( fp ) )
3414#endif
3415    {
3416        fgets( b, B-1, fp );
3417        if( !strncmp( "+==========+", b, 12 ) )
3418        {
3419            junban[count] = atoi( b+12 );
3420
3421                        pt = strchr( b, ')' ) + 1;
3422                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3423            dis[junban[count]] = (double)opt;
3424            count++;
3425
3426        }
3427                else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3428                {
3429                        break;
3430                }
3431
3432    }
3433        if( !count ) return -1;
3434
3435        count = 0;
3436    while( 1 )
3437        {
3438                if( strncmp( ">>+==========+", b, 14 ) )
3439                {
3440                        fgets( b, B-1, fp );
3441                        if( feof( fp ) ) break;
3442                        continue;
3443                }
3444                junban[count++] = atoi( b+14 );
3445//              fprintf( stderr, "t = %d\n", atoi( b+14 ) );
3446                while( fgets( b, B-1, fp ) )
3447                        if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
3448                pt = strstr( b, ":" ) +1;
3449                opt = atoi( pt );
3450
3451
3452                while( fgets( b, B-1, fp ) )
3453                        if( !strncmp( "_overlap:", b+4, 9 ) ) break;
3454                pt = strstr( b, ":" ) +1;
3455                overlapaa = atoi( pt );
3456
3457                while( fgets( b, B-1, fp ) )
3458                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3459                pt = strstr( b, ":" ) +1;
3460                qstart = atoi( pt ) - 1;
3461
3462                while( fgets( b, B-1, fp ) )
3463                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3464                pt = strstr( b, ":" ) +1;
3465                qend = atoi( pt ) - 1;
3466
3467                while( fgets( b, B-1, fp ) )
3468                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3469                pt = strstr( b, ":" ) +1;
3470                qal_display_start = atoi( pt ) - 1;
3471
3472                pt = qal;
3473                while( (c = fgetc( fp )) )
3474                {
3475                        if( c == '>' ) 
3476                        {
3477                                ungetc( c, fp );
3478                                break;
3479                        }
3480                        if( isalpha( c ) || c == '-' ) 
3481                        *pt++ = c;
3482                }
3483                *pt = 0;
3484
3485                while( fgets( b, B-1, fp ) )
3486                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3487                pt = strstr( b, ":" ) + 1;
3488                tstart = atoi( pt ) - 1;
3489
3490                while( fgets( b, B-1, fp ) )
3491                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3492                pt = strstr( b, ":" ) + 1;
3493                tend = atoi( pt ) - 1;
3494
3495                while( fgets( b, B-1, fp ) )
3496                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3497                pt = strstr( b, ":" ) + 1;
3498                tal_display_start = atoi( pt ) - 1;
3499
3500                pt = tal;
3501                while( ( c = fgetc( fp ) ) )
3502                {
3503                        if( c == '>' ) 
3504                        {
3505                                ungetc( c, fp );
3506                                break;
3507                        }
3508                        if( isalpha( c ) || c == '-' ) 
3509                        *pt++ = c;
3510                }
3511                *pt = 0;
3512
3513//              fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
3514//              fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
3515
3516//              fprintf( stderr, "qal = %s\n", qal );
3517//              fprintf( stderr, "tal = %s\n", tal );
3518
3519                qal2 = cutal( qal, qal_display_start, qstart, qend );
3520                tal2 = cutal( tal, tal_display_start, tstart, tend );
3521
3522//              fprintf( stderr, "qal2 = %s\n", qal2 );
3523//              fprintf( stderr, "tal2 = %s\n", tal2 );
3524
3525//              fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
3526                putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
3527        }
3528//      fprintf( stderr, "count = %d\n", count );
3529    return count;
3530}
3531int ReadFasta34m10_scoreonly_nucbk( FILE *fp, double *dis, int nin )
3532{
3533    int count=0;
3534    char b[B];
3535        char *pt;
3536    int pos;
3537        int opt;
3538        double z, bits;
3539
3540    count = 0;
3541    while( !feof( fp ) )
3542    {
3543        fgets( b, B-1, fp );
3544        if( !strncmp( "+===========+", b, 13 ) )
3545        {
3546            pos = atoi( b+13 );
3547
3548                        if( strchr( b, 'r' ) ) continue;
3549
3550//                      pt = strchr( b, ')' ) + 1;
3551                        pt = strchr( b, ']' ) + 1;
3552                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3553            dis[pos] += (double)opt;
3554            count++;
3555#if 0
3556                        fprintf( stderr, "b=%s\n", b );
3557                        fprintf( stderr, "opt=%d\n", opt );
3558                        fprintf( stderr, "pos=%d\n", pos );
3559                        fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3560#endif
3561
3562        }
3563                else if( 0 == strncmp( ">>><<<", b, 6 ) )
3564                {
3565                        break;
3566                }
3567
3568    }
3569        if( !count ) return -1;
3570
3571    return count;
3572}
3573
3574int ReadFasta34m10_scoreonly_nuc( FILE *fp, double *dis, int nin )
3575{
3576    int count=0;
3577    char b[B];
3578        char *pt;
3579    int pos;
3580        int opt;
3581        double z, bits;
3582        int c;
3583        int *yonda;
3584
3585
3586        yonda = AllocateIntVec( nin );
3587        for( c=0; c<nin; c++ ) yonda[c] = 0;
3588        for( c=0; c<nin; c++ ) dis[c] = 0.0;
3589
3590    count = 0;
3591    while( !feof( fp ) )
3592    {
3593        fgets( b, B-1, fp );
3594        if( !strncmp( "+===========+", b, 13 ) )
3595        {
3596            pos = atoi( b+13 );
3597
3598                        if( strchr( b, 'r' ) ) continue;
3599
3600//                      pt = strchr( b, ')' ) + 1;
3601                        pt = strchr( b, ']' ) + 1;
3602                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3603                        if( yonda[pos] == 0 )
3604                        {
3605                    dis[pos] += (double)opt;
3606                                yonda[pos] = 1;
3607                        }
3608            count++;
3609#if 0
3610                        fprintf( stderr, "b=%s\n", b );
3611                        fprintf( stderr, "opt=%d\n", opt );
3612                        fprintf( stderr, "pos=%d\n", pos );
3613                        fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3614#endif
3615
3616        }
3617        else if( !strncmp( ">>>", b, 3 ) )
3618                {
3619                        for( c=0; c<nin; c++ ) yonda[c] = 0;
3620                }
3621                else if( 0 == strncmp( ">>><<<", b, 6 ) )
3622                {
3623                        break;
3624                }
3625
3626    }
3627
3628        free( yonda );
3629
3630        if( !count ) return -1;
3631
3632    return count;
3633}
3634
3635int ReadFasta34m10_scoreonly( FILE *fp, double *dis, int nin )
3636{
3637    int count=0;
3638    char b[B];
3639        char *pt;
3640    int pos;
3641        int opt;
3642        double z, bits;
3643        int c;
3644        int *yonda;
3645
3646
3647        yonda = AllocateIntVec( nin );
3648        for( c=0; c<nin; c++ ) yonda[c] = 0;
3649        for( c=0; c<nin; c++ ) dis[c] = 0.0;
3650
3651    count = 0;
3652    while( !feof( fp ) )
3653    {
3654        fgets( b, B-1, fp );
3655        if( !strncmp( "+===========+", b, 13 ) )
3656        {
3657            pos = atoi( b+13 );
3658
3659                        pt = strchr( b, ')' ) + 1;
3660                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3661                        if( yonda[pos] == 0 )
3662                        {
3663                    dis[pos] += (double)opt;
3664                                yonda[pos] = 1;
3665                        }
3666            count++;
3667#if 0
3668                        fprintf( stderr, "b=%s\n", b );
3669                        fprintf( stderr, "opt=%d\n", opt );
3670                        fprintf( stderr, "pos=%d\n", pos );
3671                        fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3672#endif
3673
3674        }
3675        else if( !strncmp( ">>>", b, 3 ) )
3676                {
3677                        for( c=0; c<nin; c++ ) yonda[c] = 0;
3678                }
3679                else if( 0 == strncmp( ">>><<<", b, 6 ) )
3680                {
3681                        break;
3682                }
3683
3684    }
3685
3686        free( yonda );
3687
3688        if( !count ) return -1;
3689
3690    return count;
3691}
3692int ReadFasta34( FILE *fp, double *dis, int nseq, char name[M][B], LocalHom *localhomlist )
3693{
3694    int count=0;
3695    char b[B];
3696        char *pt;
3697    static int junban[M];
3698        int overlapaa;
3699        int opt, qstart, qend, tstart, tend;
3700        double z, bits;
3701
3702
3703    count = 0;
3704#if 0
3705    for( i=0; i<10000000 && count<nseq; i++ )
3706#else
3707    while( !feof( fp ) )
3708#endif
3709    {
3710        fgets( b, B-1, fp );
3711        if( !strncmp( "+==========+", b, 12 ) )
3712        {
3713            junban[count] = atoi( b+12 );
3714
3715                        pt = strchr( b, ')' ) + 1;
3716                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3717            dis[junban[count]] = (double)opt;
3718            count++;
3719
3720        }
3721                else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3722                {
3723                        break;
3724                }
3725
3726    }
3727        if( !count ) return -1;
3728
3729        count = 0;
3730    while( !feof( fp ) )
3731        {
3732                if( !strncmp(">>+==========+", b, 14 ) )
3733                {
3734            junban[count] = atoi( b+14 );
3735            count++;
3736                fgets( b, B-1, fp ); // initn:
3737                        pt = strstr( b, "opt: " ) + 5;
3738                        localhomlist[junban[count-1]].opt = atof( pt );
3739                fgets( b, B-1, fp ); // Smith-Waterman score
3740                        pt = strstr( b, "ungapped) in " ) + 13;
3741                        sscanf( pt, "%d", &overlapaa ); 
3742                        fprintf( stderr, "pt = %s, overlapaa = %d\n", pt, overlapaa );
3743                        pt = strstr( b, "overlap (" ) + 8;
3744                        sscanf( pt, "(%d-%d:%d-%d)", &qstart, &qend, &tstart, &tend ); 
3745                        localhomlist[junban[count-1]].overlapaa = overlapaa;
3746                        localhomlist[junban[count-1]].start1 = qstart-1;
3747                        localhomlist[junban[count-1]].end1   = qend-1;
3748                        localhomlist[junban[count-1]].start2 = tstart-1;
3749                        localhomlist[junban[count-1]].end2   = tend-1;
3750                }
3751        fgets( b, B-1, fp );
3752        }
3753        fprintf( stderr, "count = %d\n", count );
3754    return count;
3755}
3756
3757int ReadFasta3( FILE *fp, double *dis, int nseq, char name[M][B] )
3758{
3759    int count=0;
3760    char b[B];
3761        char *pt;
3762    int junban[M];
3763        int initn, init1, opt;
3764        double z;
3765
3766    count = 0;
3767#if 0
3768    for( i=0; i<10000000 && count<nseq; i++ )
3769#else
3770    while( !feof( fp ) )
3771#endif
3772    {
3773        fgets( b, B-1, fp );
3774        if( !strncmp( "+==========+", b, 12 ) )
3775        {
3776            junban[count] = atoi( b+12 );
3777
3778                        pt = strchr( b, ')' ) + 1;
3779                        sscanf( pt, "%d %d %d %lf", &initn, &init1, &opt, &z ); 
3780            dis[junban[count]] = (double)opt;
3781            count++;
3782        }
3783    }
3784    return 0;
3785}
3786
3787int ReadFasta( FILE *fp, double *dis, int nseq, char name[M][B] )
3788{
3789    int i, count=0;
3790    char b[B];
3791    int junban[M];
3792        int initn, init1, opt;
3793
3794    count = 0;
3795        for( i=0; i<nseq; i++ ) dis[i] = 0.0;
3796    for( i=0; !feof( fp ) && count<nseq; i++ )
3797    {
3798        fgets( b, B-1, fp );
3799        if( !strncmp( "+==========+", b, 12 ) )
3800        {
3801            junban[count] = atoi( b+12 );
3802
3803                        sscanf( b+50, "%d %d %d", &initn, &init1, &opt ); 
3804            dis[junban[count]] = (double)opt;
3805            count++;
3806        }
3807    }
3808
3809/*
3810    count = 0;
3811    for( i=0; i<100000 && count<nseq; i++ )
3812    {
3813        fgets( b, B-1, fp );
3814        if( !strncmp( name[junban[count]], b, 20  ) )
3815        {
3816            dis[junban[count]] = atof( b+65 );
3817            count++;
3818        }
3819    }
3820*/
3821    return 0;
3822}
3823
3824
3825int ReadOpt( FILE *fp, int opt[M], int nseq, char name[M][B] )
3826{
3827    int i, count=0;
3828    char b[B];
3829    int junban[M];
3830        int optt, initn, init1;
3831
3832    count = 0;
3833    for( i=0; i<10000000 && count<nseq; i++ )
3834    {
3835        fgets( b, B-1, fp );
3836        if( !strncmp( "+==========+", b, 12 ) )
3837        {
3838            junban[count] = atoi( b+12 );
3839                        sscanf( b+50, "%d %d %d", &initn, &init1, &optt ); 
3840            opt[junban[count]] = (double)optt;
3841            count++;
3842        }
3843    }
3844    return 0;
3845}
3846
3847int ReadOpt2( FILE *fp, int opt[M], int nseq, char name[M][B] )
3848{
3849    int i, count=0;
3850    char b[B];
3851    int junban[M];
3852
3853    count = 0;
3854    for( i=0; i<10000000 && count<nseq; i++ )
3855    {
3856        fgets( b, B-1, fp );
3857        if( !strncmp( "+==========+", b, 12 ) )
3858        {
3859            junban[count] = atoi( b+12 );
3860            opt[junban[count]] = atoi( b+65 );
3861            count++;
3862        }
3863    }
3864    return 0;
3865}
3866
3867
3868
3869int writePre( int nseq, char **name, int nlen[M], char **aseq, int force )
3870{
3871#if USE_XCED
3872        int i, value;
3873        if( !signalSM )
3874        {
3875                if( force ) 
3876                {
3877                        rewind( prep_g );
3878                        if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3879#if 0
3880                        else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3881#else
3882                        else    writeData( prep_g, nseq, name, nlen, aseq );
3883#endif
3884                }
3885                return( 0 );
3886        }
3887        for( i=0; i<10; i++ )
3888        {
3889#if IODEBUG
3890                fprintf( stderr, "SEMAPHORE = %d\n", signalSM[SEMAPHORE] );
3891#endif
3892                if( signalSM[SEMAPHORE]-- > 0 )
3893                {
3894#if 0 /* /tmp/pre €ÎŽØ·ž€Ç€Ï€º€·€¿ */
3895                        if( ferror( prep_g ) ) prep_g = fopen( "pre", "w" );
3896                        if( !prep_g ) ErrorExit( "Cannot re-open pre." );
3897#endif
3898                        rewind( prep_g );
3899                        signalSM[STATUS] = IMA_KAITERU;
3900#if IODEBUG
3901                        if( force ) fprintf( stderr, "FINAL " );
3902#endif
3903                        if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3904                        else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3905                        /*
3906                        fprintf( prep_g, '\EOF' );
3907                        */
3908                        fflush( prep_g );
3909                        if( force ) signalSM[STATUS] = OSHIMAI;
3910                        else        signalSM[STATUS] = KAKIOWATTA;
3911                        value = 1;
3912                        signalSM[SEMAPHORE]++;
3913#if IODEBUG
3914                        fprintf( stderr, "signalSM[STATUS] = %c\n", signalSM[STATUS] );
3915#endif
3916                        break;
3917                }
3918                else
3919                {
3920#if IODEBUG
3921                        fprintf( stderr, "YONDERUKARA_AKIRAMERU\n" );
3922#endif
3923                        value = 0;
3924                        signalSM[SEMAPHORE]++;
3925                        if( !force ) break;
3926#if IODEBUG
3927                        fprintf( stderr, "MATSU\n" );
3928#endif
3929                        sleep( 1 );
3930                }
3931        }
3932        if( force && !value ) ErrorExit( "xced ga pre wo hanasanai \n" );
3933        return( value );
3934#else
3935        if( force ) 
3936        {
3937                rewind( prep_g );
3938                writeData_pointer( prep_g, nseq, name, nlen, aseq );
3939        }
3940#endif
3941        return( 0 );
3942}
3943
3944
3945void readOtherOptions( int *ppidptr, int *fftThresholdptr, int *fftWinSizeptr )
3946{
3947        if( calledByXced )
3948        {
3949                FILE *fp = fopen( "pre", "r" );
3950                char b[B];
3951                if( !fp ) ErrorExit( "Cannot open pre.\n" );
3952                fgets( b, B-1, fp );
3953                sscanf( b, "%d %d %d", ppidptr, fftThresholdptr, fftWinSizeptr );
3954                fclose( fp );
3955#if IODEBUG
3956        fprintf( stderr, "b = %s\n", b );
3957        fprintf( stderr, "ppid = %d\n", ppid );
3958        fprintf( stderr, "fftThreshold = %d\n", fftThreshold );
3959        fprintf( stderr, "fftWinSize = %d\n", fftWinSize );
3960#endif
3961        }
3962        else
3963        {
3964                *ppidptr = 0;
3965                *fftThresholdptr = FFT_THRESHOLD;
3966                if( dorp == 'd' )
3967                        *fftWinSizeptr = FFT_WINSIZE_D;
3968                else
3969                        *fftWinSizeptr = FFT_WINSIZE_P;
3970        }
3971#if 0
3972        fprintf( stderr, "fftThresholdptr=%d\n", *fftThresholdptr );
3973        fprintf( stderr, "fftWinSizeptr=%d\n", *fftWinSizeptr );
3974#endif
3975}
3976
3977void initSignalSM( void )
3978{
3979//      int signalsmid;
3980
3981#if IODEBUG
3982        if( ppid ) fprintf( stderr, "PID of xced = %d\n", ppid );
3983#endif
3984        if( !ppid )
3985        {
3986                signalSM = NULL;
3987                return;
3988        }
3989
3990#if 0
3991        signalsmid = shmget( (key_t)ppid, 3, IPC_ALLOC | 0666 );
3992        if( signalsmid == -1 ) ErrorExit( "Cannot get Shared memory for signal.\n" );
3993        signalSM = shmat( signalsmid, 0, 0 );
3994        if( (int)signalSM == -1 ) ErrorExit( "Cannot attatch Shared Memory for signal!\n" );
3995        signalSM[STATUS] = IMA_KAITERU;
3996        signalSM[SEMAPHORE] = 1;
3997#endif
3998}
3999
4000void initFiles( void )
4001{
4002        char pname[100];
4003        if( ppid )
4004                sprintf( pname, "/tmp/pre.%d", ppid );
4005        else
4006                sprintf( pname, "pre" );
4007        prep_g = fopen( pname, "w" );
4008        if( !prep_g ) ErrorExit( "Cannot open pre" );
4009
4010        trap_g = fopen( "trace", "w" );
4011        if( !trap_g ) ErrorExit( "cannot open trace" );
4012        fprintf( trap_g, "PID = %d\n", getpid() );
4013        fflush( trap_g );
4014}
4015
4016
4017void WriteForFasta( FILE *fp, int locnjob, char **name, int nlen[M], char **aseq )
4018{
4019    static char b[N];
4020    int i, j;
4021    int nalen[M];
4022
4023    for( i=0; i<locnjob; i++ )
4024    {
4025        nalen[i] = strlen( aseq[i] );
4026        fprintf( fp, ">%s\n", name[i] );
4027        for( j=0; j<nalen[i]; j=j+C ) 
4028        {
4029            strncpy( b, aseq[i]+j, C ); b[C] = 0;
4030            fprintf( fp, "%s\n",b );
4031        }
4032    }
4033}
4034
4035void readlocalhomtable2( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
4036{
4037        double opt;
4038        static char buff[B];
4039        char infor[100];
4040        int i, j, overlapaa, start1, end1, start2, end2;
4041        LocalHom *tmpptr1, *tmpptr2;
4042
4043//      for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
4044
4045        while ( NULL != fgets( buff, B-1, fp ) )
4046        {
4047//              fprintf( stderr, "\n" );
4048                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4049                if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
4050
4051#if 0
4052                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4053#endif
4054
4055//              if( i < j )
4056                {
4057                        if( localhomtable[i][j].nokori++ > 0 )
4058                        {
4059                                tmpptr1 = localhomtable[i][j].last;
4060//                              fprintf( stderr, "reallocating, localhomtable[%d][%d].nokori = %d\n", i, j, localhomtable[i][j].nokori );
4061                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4062                                tmpptr1 = tmpptr1->next;
4063                                tmpptr1->extended = -1;
4064                                tmpptr1->next = NULL;
4065                                localhomtable[i][j].last = tmpptr1;
4066//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
4067                        }
4068                        else
4069                        {
4070                                tmpptr1 = localhomtable[i]+j;
4071//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
4072                        }
4073       
4074                        tmpptr1->start1 = start1;
4075                        tmpptr1->start2 = start2;
4076                        tmpptr1->end1 = end1;
4077                        tmpptr1->end2 = end2;
4078//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4079//                      tmpptr1->opt = opt;
4080                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4081                        tmpptr1->overlapaa = overlapaa;
4082                        tmpptr1->korh = *infor;
4083                }
4084//              else
4085                {
4086                        if( localhomtable[j][i].nokori++ > 0 )
4087                        {
4088                                tmpptr2 = localhomtable[j][i].last;
4089                                tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4090                                tmpptr2 = tmpptr2->next;
4091                                tmpptr2->extended = -1;
4092                                tmpptr2->next = NULL;
4093                                localhomtable[j][i].last = tmpptr2;
4094//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
4095                        }
4096                        else
4097                        {
4098                                tmpptr2 = localhomtable[j]+i;
4099//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
4100                        }
4101       
4102                        tmpptr2->start2 = start1;
4103                        tmpptr2->start1 = start2;
4104                        tmpptr2->end2 = end1;
4105                        tmpptr2->end1 = end2;
4106//                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4107//                      tmpptr2->opt = opt;
4108                        tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4109                        tmpptr2->overlapaa = overlapaa;
4110                        tmpptr2->korh = *infor;
4111       
4112//                      fprintf( stderr, "i=%d, j=%d, st1=%d, en1=%d, opt = %f\n", i, j, tmpptr1->start1, tmpptr1->end1, opt );
4113                }
4114
4115        }
4116}
4117void readlocalhomtable( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
4118{
4119        double opt;
4120        static char buff[B];
4121        char infor[100];
4122        int i, j, overlapaa, start1, end1, start2, end2;
4123        int **nlocalhom = NULL;
4124        LocalHom *tmpptr1=NULL, *tmpptr2=NULL; // by D.Mathog, a guess
4125
4126        nlocalhom = AllocateIntMtx( njob, njob );
4127        for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
4128
4129        while ( NULL != fgets( buff, B-1, fp ) )
4130        {
4131//              fprintf( stderr, "\n" );
4132                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4133                if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
4134
4135#if 0
4136                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4137#endif
4138
4139
4140//              if( i < j )
4141                {
4142                        if( nlocalhom[i][j]++ > 0 )
4143                        {
4144//                              fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4145                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4146                                tmpptr1 = tmpptr1->next;
4147                                tmpptr1->next = NULL;
4148                        }
4149                        else
4150                        {
4151                                tmpptr1 = localhomtable[i]+j;
4152//                              fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4153                        }
4154       
4155                        tmpptr1->start1 = start1; // CHUUI!!!!
4156                        tmpptr1->start2 = start2;
4157                        tmpptr1->end1 = end1; // CHUUI!!!!
4158                        tmpptr1->end2 = end2;
4159//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4160//                      tmpptr1->opt = opt;
4161                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4162                        tmpptr1->overlapaa = overlapaa;
4163                        tmpptr1->korh = *infor;
4164       
4165//                      fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4166                }
4167//              else
4168                {
4169                        if( nlocalhom[j][i]++ > 0 )
4170                        {
4171                                tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4172                                tmpptr2 = tmpptr2->next;
4173                                tmpptr2->next = NULL;
4174                        }
4175                        else
4176                                tmpptr2 = localhomtable[j]+i;
4177       
4178                        tmpptr2->start2 = start1; // CHUUI!!!!
4179                        tmpptr2->start1 = start2;
4180                        tmpptr2->end2 = end1; // CHUUI!!!!
4181                        tmpptr2->end1 = end2;
4182//                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4183//                      tmpptr2->opt = opt;
4184                        tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4185                        tmpptr2->overlapaa = overlapaa;
4186                        tmpptr2->korh = *infor;
4187
4188//                      fprintf( stderr, "j=%d, i=%d, opt = %f\n", j, i, opt );
4189                }
4190
4191        }
4192        FreeIntMtx( nlocalhom );
4193}
4194
4195
4196void readlocalhomtable_two( FILE*fp, int norg, int nadd, LocalHom **localhomtable, LocalHom **localhomtablex, char *kozoarivec ) // for test only
4197{
4198        double opt;
4199        static char buff[B];
4200        char infor[100];
4201        int i, j, overlapaa, start1, end1, start2, end2;
4202        int **nlocalhom = NULL;
4203        int **nlocalhomx = NULL;
4204        LocalHom *tmpptr1=NULL, *tmpptr2=NULL; // by D.Mathog, a guess
4205
4206        nlocalhom = AllocateIntMtx( norg, nadd );
4207        for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ ) nlocalhom[i][j] = 0;
4208        nlocalhomx = AllocateIntMtx( nadd, norg );
4209        for( i=0; i<nadd; i++ ) for( j=0; j<norg; j++ ) nlocalhomx[i][j] = 0;
4210
4211        while ( NULL != fgets( buff, B-1, fp ) )
4212        {
4213//              fprintf( stderr, "\n" );
4214                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4215                if( *infor == 'k' ) 
4216                {
4217                        fprintf( stderr, "Not supported!\n" );
4218                        exit( 1 );
4219                }
4220                j -= norg;
4221
4222#if 0
4223                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4224#endif
4225
4226
4227//              if( i < j )
4228                {
4229                        if( nlocalhom[i][j]++ > 0 )
4230                        {
4231//                              fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4232                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4233                                tmpptr1 = tmpptr1->next;
4234                                tmpptr1->next = NULL;
4235                        }
4236                        else
4237                        {
4238                                tmpptr1 = localhomtable[i]+j;
4239//                              fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4240                        }
4241       
4242                        tmpptr1->start1 = start1; // CHUUI!!!!
4243                        tmpptr1->start2 = start2;
4244                        tmpptr1->end1 = end1; // CHUUI!!!!
4245                        tmpptr1->end2 = end2;
4246//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4247//                      tmpptr1->opt = opt;
4248                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4249                        tmpptr1->overlapaa = overlapaa;
4250                        tmpptr1->korh = *infor;
4251       
4252//                      fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4253                }
4254
4255                {
4256                        if( nlocalhomx[j][i]++ > 0 )
4257                        {
4258                                tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4259                                tmpptr2 = tmpptr2->next;
4260                                tmpptr2->next = NULL;
4261                        }
4262                        else
4263                                tmpptr2 = localhomtablex[j]+i;
4264       
4265                        tmpptr2->start2 = start1+1; // CHUUI!!!!
4266                        tmpptr2->start1 = start2;
4267                        tmpptr2->end2 = end1+1; // CHUUI!!!!
4268                        tmpptr2->end1 = end2;
4269//                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4270//                      tmpptr2->opt = opt;
4271                        tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4272                        tmpptr2->overlapaa = overlapaa;
4273                        tmpptr2->korh = *infor;
4274
4275//                      fprintf( stderr, "j=%d, i=%d, opt = %f\n", j, i, opt );
4276                }
4277
4278        }
4279        FreeIntMtx( nlocalhom );
4280        FreeIntMtx( nlocalhomx );
4281}
4282
4283void readlocalhomtable_one( FILE*fp, int norg, int nadd, LocalHom **localhomtable, char *kozoarivec ) // for test only
4284{
4285        double opt;
4286        static char buff[B];
4287        char infor[100];
4288        int i, j, overlapaa, start1, end1, start2, end2;
4289        int **nlocalhom = NULL;
4290        LocalHom *tmpptr1=NULL; // by D.Mathog, a guess
4291
4292        nlocalhom = AllocateIntMtx( norg, nadd );
4293        for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ ) nlocalhom[i][j] = 0;
4294
4295        while ( NULL != fgets( buff, B-1, fp ) )
4296        {
4297//              fprintf( stderr, "\n" );
4298                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4299                if( *infor == 'k' ) 
4300                {
4301                        fprintf( stderr, "Not supported!\n" );
4302                        exit( 1 );
4303                }
4304                j -= norg;
4305
4306#if 0
4307                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4308#endif
4309
4310
4311//              if( i < j )
4312                {
4313                        if( nlocalhom[i][j]++ > 0 )
4314                        {
4315//                              fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4316                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4317                                tmpptr1 = tmpptr1->next;
4318                                tmpptr1->next = NULL;
4319                        }
4320                        else
4321                        {
4322                                tmpptr1 = localhomtable[i]+j;
4323//                              fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4324                        }
4325       
4326                        tmpptr1->start1 = start1; // CHUUI!!!!
4327                        tmpptr1->start2 = start2;
4328                        tmpptr1->end1 = end1; // CHUUI!!!!
4329                        tmpptr1->end2 = end2;
4330//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4331//                      tmpptr1->opt = opt;
4332                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4333                        tmpptr1->overlapaa = overlapaa;
4334                        tmpptr1->korh = *infor;
4335       
4336//                      fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4337                }
4338
4339        }
4340        FreeIntMtx( nlocalhom );
4341}
4342
4343void outlocalhom_part( LocalHom **localhom, int norg, int nadd )
4344{
4345        int i, j;
4346        LocalHom *tmpptr;
4347        for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ )
4348        {
4349                tmpptr = localhom[i]+j;
4350                fprintf( stderr, "%d-%d\n", i, j+norg );
4351                do
4352                {
4353                        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 );
4354                }
4355                while( (tmpptr=tmpptr->next) );
4356        }
4357}
4358
4359void outlocalhom( LocalHom **localhom, int nseq )
4360{
4361        int i, j;
4362        LocalHom *tmpptr;
4363        for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
4364        {
4365                tmpptr = localhom[i]+j;
4366                fprintf( stderr, "%d-%d\n", i, j );
4367                do
4368                {
4369                        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 );
4370                }
4371                while( (tmpptr=tmpptr->next) );
4372        }
4373}
4374
4375void outlocalhompt( LocalHom ***localhom, int n1, int n2 )
4376{
4377        int i, j;
4378        LocalHom *tmpptr;
4379        for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
4380        {
4381                tmpptr = localhom[i][j];
4382//              fprintf( stdout, "%d-%d\n", i, j );
4383                do
4384                {
4385                        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 );
4386                }
4387                while( (tmpptr=tmpptr->next) );
4388        }
4389}
4390
4391void FreeLocalHomTable_part( LocalHom **localhomtable, int n, int m ) 
4392{
4393        int i, j;
4394        LocalHom *ppp, *tmpptr;
4395        for( i=0; i<n; i++ ) 
4396        {
4397                for( j=0; j<m; j++ )
4398                {
4399                        tmpptr=localhomtable[i]+j;
4400                        ppp = tmpptr->next;
4401                        for( ; tmpptr; tmpptr=ppp )
4402                        {
4403#if DEBUG
4404                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4405#endif
4406                                ppp = tmpptr->next;
4407                                if( tmpptr!=localhomtable[i]+j ) 
4408                                {
4409#if DEBUG
4410                                        fprintf( stderr, "freeing %p\n", tmpptr );
4411#endif
4412                                        free( tmpptr );
4413                                }
4414                        }
4415                }
4416#if DEBUG
4417                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4418#endif
4419                free( localhomtable[i] );
4420        }
4421#if DEBUG
4422        fprintf( stderr, "freeing localhomtable\n" );
4423#endif
4424        free( localhomtable );
4425#if DEBUG
4426        fprintf( stderr, "freed\n" );
4427#endif
4428}
4429
4430void FreeLocalHomTable_two( LocalHom **localhomtable, int n, int m ) 
4431{
4432        int i, j;
4433        LocalHom *ppp, *tmpptr;
4434        for( i=0; i<n; i++ ) 
4435        {
4436                for( j=0; j<m; j++ )
4437                {
4438                        tmpptr=localhomtable[i]+j;
4439                        ppp = tmpptr->next;
4440                        for( ; tmpptr; tmpptr=ppp )
4441                        {
4442#if DEBUG
4443                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4444#endif
4445                                ppp = tmpptr->next;
4446                                if( tmpptr!=localhomtable[i]+j ) 
4447                                {
4448#if DEBUG
4449                                        fprintf( stderr, "freeing %p\n", tmpptr );
4450#endif
4451                                        free( tmpptr );
4452                                }
4453                        }
4454                }
4455#if DEBUG
4456                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4457#endif
4458                free( localhomtable[i] );
4459        }
4460
4461        for( i=n; i<n+m; i++ ) 
4462        {
4463                for( j=0; j<n; j++ )
4464                {
4465                        tmpptr=localhomtable[i]+j;
4466                        ppp = tmpptr->next;
4467                        for( ; tmpptr; tmpptr=ppp )
4468                        {
4469#if DEBUG
4470                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4471#endif
4472                                ppp = tmpptr->next;
4473                                if( tmpptr!=localhomtable[i]+j ) 
4474                                {
4475#if DEBUG
4476                                        fprintf( stderr, "freeing %p\n", tmpptr );
4477#endif
4478                                        free( tmpptr );
4479                                }
4480                        }
4481                }
4482#if DEBUG
4483                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4484#endif
4485                free( localhomtable[i] );
4486        }
4487#if DEBUG
4488        fprintf( stderr, "freeing localhomtable\n" );
4489#endif
4490        free( localhomtable );
4491#if DEBUG
4492        fprintf( stderr, "freed\n" );
4493#endif
4494}
4495
4496void FreeLocalHomTable_one( LocalHom **localhomtable, int n, int m ) 
4497{
4498        int i, j;
4499        LocalHom *ppp, *tmpptr;
4500        for( i=0; i<n; i++ ) 
4501        {
4502                for( j=0; j<m; j++ )
4503                {
4504                        tmpptr=localhomtable[i]+j;
4505                        ppp = tmpptr->next;
4506                        for( ; tmpptr; tmpptr=ppp )
4507                        {
4508#if DEBUG
4509                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4510#endif
4511                                ppp = tmpptr->next;
4512                                if( tmpptr!=localhomtable[i]+j ) 
4513                                {
4514#if DEBUG
4515                                        fprintf( stderr, "freeing %p\n", tmpptr );
4516#endif
4517                                        free( tmpptr );
4518                                }
4519                        }
4520                }
4521#if DEBUG
4522                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4523#endif
4524                free( localhomtable[i] );
4525        }
4526
4527#if DEBUG
4528        fprintf( stderr, "freeing localhomtable\n" );
4529#endif
4530        free( localhomtable );
4531#if DEBUG
4532        fprintf( stderr, "freed\n" );
4533#endif
4534}
4535
4536void FreeLocalHomTable( LocalHom **localhomtable, int n ) 
4537{
4538        int i, j;
4539        LocalHom *ppp, *tmpptr;
4540        for( i=0; i<n; i++ ) 
4541        {
4542                for( j=0; j<n; j++ )
4543                {
4544                        tmpptr=localhomtable[i]+j;
4545                        ppp = tmpptr->next;
4546                        for( ; tmpptr; tmpptr=ppp )
4547                        {
4548#if DEBUG
4549                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4550#endif
4551                                ppp = tmpptr->next;
4552                                if( tmpptr!=localhomtable[i]+j ) 
4553                                {
4554#if DEBUG
4555                                        fprintf( stderr, "freeing %p\n", tmpptr );
4556#endif
4557                                        free( tmpptr );
4558                                }
4559                        }
4560                }
4561#if DEBUG
4562                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4563#endif
4564                free( localhomtable[i] );
4565        }
4566#if DEBUG
4567        fprintf( stderr, "freeing localhomtable\n" );
4568#endif
4569        free( localhomtable );
4570#if DEBUG
4571        fprintf( stderr, "freed\n" );
4572#endif
4573}
4574
4575char *progName( char *str )
4576{
4577    char *value; 
4578    if( ( value = strrchr( str, '/' ) ) != NULL )
4579        return( value+1 );
4580    else   
4581        return( str );
4582}
4583
4584static void tabtospace( char *str )
4585{
4586        char *p;
4587//      fprintf( stderr, "before = %s\n", str );
4588        while( NULL != ( p = strchr( str , '\t' ) ) )
4589        {
4590                *p = ' ';
4591        }
4592//      fprintf( stderr, "after = %s\n", str );
4593}
4594
4595static char *extractfirstword( char *str )
4596{
4597        char *val = str;
4598
4599        tabtospace( str );
4600        while( *str )
4601        {
4602                if( val == str && *str == ' ' )
4603                {
4604                        val++; str++;
4605                }
4606                else if( *str != ' ' )
4607                {
4608                        str++;
4609                }
4610                else if( *str == ' ' )
4611                {
4612                        *str = 0;
4613                }
4614        }
4615        return( val );
4616}
4617
4618void phylipout_pointer( FILE *fp, int nseq, int maxlen, char **seq, char **name, int *order, int namelen )
4619{
4620        int pos, pos2, j;
4621        if( namelen == -1 ) namelen = 10;
4622        pos = 0;
4623
4624        fprintf( fp, " %d %d\n", nseq, maxlen );
4625       
4626        while( pos < maxlen )
4627        {
4628                for( j=0; j<nseq; j++ )
4629                {
4630                        if( pos == 0 )
4631                                fprintf( fp, "%-*.*s", namelen, namelen, extractfirstword( name[order[j]]+1 ) );
4632                        else
4633                                fprintf( fp, "%-*.*s", namelen, namelen, "" );
4634
4635                        pos2 = pos;
4636                        while( pos2 < pos+41 && pos2 < maxlen )
4637                        {
4638                                fprintf( fp, " %.10s", seq[order[j]]+pos2 );
4639                                pos2 += 10;
4640                        }
4641                        fprintf( fp, "\n" );
4642                }
4643                fprintf( fp, "\n" );
4644                pos += 50;
4645        }
4646}
4647
4648void clustalout_pointer( FILE *fp, int nseq, int maxlen, char **seq, char **name, char *mark, char *comment, int *order, int namelen )
4649{
4650        int pos, j;
4651        if( namelen == -1 ) namelen = 15;
4652        pos = 0;
4653        if( comment == NULL )
4654                fprintf( fp, "CLUSTAL format alignment by MAFFT (v%s)\n\n", VERSION );
4655        else
4656                fprintf( fp, "CLUSTAL format alignment by MAFFT %s (v%s)\n\n", comment, VERSION );
4657       
4658        while( pos < maxlen )
4659        {
4660                fprintf( fp, "\n" );
4661                for( j=0; j<nseq; j++ )
4662                {
4663                        fprintf( fp, "%-*.*s ", namelen, namelen, extractfirstword( name[order[j]]+1 ) );
4664                        fprintf( fp, "%.60s\n", seq[order[j]]+pos ); // Ĺ€µ€¬°ã€Š€È€À€á¡£
4665                }
4666                if( mark )
4667                {
4668                        fprintf( fp, "%-*.*s ", namelen, namelen, "" );
4669                        fprintf( fp, "%.60s\n", mark + pos ); // Ĺ€µ€¬°ã€Š€È€À€á¡£
4670                }
4671                pos += 60;
4672        }
4673}
4674
4675
4676void writeData_reorder_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq, int *order )
4677{
4678        int i, j, k;
4679        int nalen;
4680
4681        for( i=0; i<locnjob; i++ )
4682        {
4683                k = order[i];
4684#if DEBUG
4685                fprintf( stderr, "i = %d in writeData\n", i );
4686#endif
4687                nalen = strlen( aseq[k] );
4688                fprintf( fp, ">%s\n", name[k]+1 );
4689                for( j=0; j<nalen; j=j+C )
4690                {
4691#if 0
4692                        strncpy( b, aseq[k]+j, C ); b[C] = 0;
4693                        fprintf( fp, "%s\n",b );
4694#else
4695                        fprintf( fp, "%.*s\n", C, aseq[k]+j );
4696#endif
4697                }
4698        }
4699}
4700void writeData_reorder( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq, int *order )
4701{
4702        int i, j, k;
4703        int nalen;
4704
4705        for( i=0; i<locnjob; i++ )
4706        {
4707                k = order[i];
4708#if DEBUG
4709                fprintf( stderr, "i = %d in writeData\n", i );
4710#endif
4711                nalen = strlen( aseq[k] );
4712                fprintf( fp, ">%s\n", name[k]+1 );
4713                for( j=0; j<nalen; j=j+C )
4714                {
4715#if 0
4716                        strncpy( b, aseq[k]+j, C ); b[C] = 0;
4717                        fprintf( fp, "%s\n",b );
4718#else
4719                        fprintf( fp, "%.*s\n", C, aseq[k]+j );
4720#endif
4721                }
4722        }
4723}
4724static void showaamtxexample()
4725{
4726        fprintf( stderr, "Format error in aa matrix\n" );
4727        fprintf( stderr, "# Example:\n" );
4728        fprintf( stderr, "# comment\n" );
4729        fprintf( stderr, "   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V\n" );
4730        fprintf( stderr, "A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0\n" );
4731        fprintf( stderr, "R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3\n" );
4732        fprintf( stderr, "...\n" );
4733        fprintf( stderr, "V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4\n" );
4734        fprintf( stderr, "frequency 0.07 0.05 0.04 0.05 0.02 .. \n" );
4735        fprintf( stderr, "# Example end\n" );
4736        fprintf( stderr, "Only the lower half is loaded\n" );
4737        fprintf( stderr, "The last line (frequency) is optional.\n" );
4738        exit( 1 );
4739}
4740
4741double *loadaamtx( void )
4742{
4743        int i, j, k, ii, jj;
4744        double *val;
4745        double **raw;
4746        int *map;
4747        char *aaorder = "ARNDCQEGHILKMFPSTWYV";
4748        char *inorder;
4749        char *line;
4750        char *ptr1;
4751        char *ptr2;
4752        char *mtxfname = "_aamtx";
4753        FILE *mf;
4754
4755        raw = AllocateDoubleMtx( 21, 20 );
4756        val = AllocateDoubleVec( 420 );
4757        map = AllocateIntVec( 20 );
4758
4759        if( dorp != 'p' )
4760        {
4761                fprintf( stderr, "User-defined matrix is not supported for DNA\n" );
4762                exit( 1 );
4763        }
4764
4765        mf = fopen( mtxfname, "r" );
4766        if( mf == NULL )
4767        {
4768                fprintf( stderr, "Cannot open the _aamtx file\n" );
4769                exit( 1 );
4770        }
4771
4772        inorder = calloc( 1000, sizeof( char ) );
4773        line = calloc( 1000, sizeof( char ) );
4774       
4775
4776        while( !feof( mf ) )
4777        {
4778                fgets( inorder, 999, mf );
4779                if( inorder[0] != '#' ) break;
4780        }
4781        ptr1 = ptr2 = inorder;
4782        while( *ptr2 )
4783        {
4784                if( isalpha( *ptr2 ) )
4785                {
4786                        *ptr1 = toupper( *ptr2 );
4787                        ptr1++;
4788                }
4789                ptr2++;
4790        }
4791        inorder[20] = 0;
4792
4793        for( i=0; i<20; i++ )
4794        {
4795                ptr2 = strchr( inorder, aaorder[i] );
4796                if( ptr2 == NULL )
4797                {
4798                        fprintf( stderr, "%c: not found in the first 20 letters.\n", aaorder[i] );
4799                        showaamtxexample();
4800                }
4801                else
4802                {
4803                        map[i] = ptr2 - inorder;
4804                }
4805        }
4806
4807        i = 0;
4808        while( !feof( mf ) )
4809        {
4810                fgets( line, 999, mf );
4811//              fprintf( stderr, "line = %s\n", line );
4812                if( line[0] == '#' ) continue;
4813                ptr1 = line;
4814//              fprintf( stderr, "line = %s\n", line );
4815                for( j=0; j<=i; j++ )
4816                {
4817                        while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4818                                ptr1++;
4819
4820                        raw[i][j] = atof( ptr1 );
4821//                      fprintf( stderr, "raw[][]=%f, %c-%c %d-%d\n", raw[i][j], inorder[i], inorder[j], i, j );
4822                        ptr1 = strchr( ptr1, ' ' );
4823                        if( ptr1 == NULL && j<i) showaamtxexample();
4824                }
4825                i++;
4826                if( i > 19 ) break;
4827        }
4828
4829        for( i=0; i<20; i++ ) raw[20][i] = -1.0;
4830        while( !feof( mf ) )
4831        {
4832                fgets( line, 999, mf );
4833                if( line[0] == 'f' )
4834                {
4835//                      fprintf( stderr, "line = %s\n", line );
4836                        ptr1 = line;
4837                        for( j=0; j<20; j++ )
4838                        {
4839                                while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4840                                        ptr1++;
4841       
4842                                raw[20][j] = atof( ptr1 );
4843//                              fprintf( stderr, "raw[20][]=%f, %c %d\n", raw[20][j], inorder[i], j );
4844                                ptr1 = strchr( ptr1, ' ' );
4845                                if( ptr1 == NULL && j<19) showaamtxexample();
4846                        }
4847                        break;
4848                }
4849        }
4850
4851        k = 0;
4852        for( i=0; i<20; i++ )
4853        {
4854                for( j=0; j<=i; j++ )
4855                {
4856                        if( i != j )
4857                        {
4858                                ii = MAX( map[i], map[j] );
4859                                jj = MIN( map[i], map[j] );
4860                        }
4861                        else ii = jj = map[i];
4862                        val[k++] = raw[ii][jj];
4863//                      fprintf( stderr, "%c-%c, %f\n", aaorder[i], aaorder[j], val[k-1] );
4864                }
4865        }
4866        for( i=0; i<20; i++ ) val[400+i] = raw[20][map[i]];
4867
4868        fprintf( stderr, "inorder = %s\n", inorder );
4869        fclose( mf );
4870        free( inorder );
4871        free( line );
4872        FreeDoubleMtx( raw );
4873        free( map );
4874        return( val );
4875}
4876
4877static void tab2space( char *s ) // nen no tame
4878{
4879        while( *s )
4880        {
4881                if( *s == '\t' ) *s = ' ';
4882                s++;
4883        }
4884}
4885
4886static int readasubalignment( char *s, int *t, int *preservegaps )
4887{
4888        int v = 0;
4889        char status = 's';
4890        char *pt = s;
4891        *preservegaps = 0;
4892        tab2space( s );
4893        while( *pt )
4894        {
4895                if( *pt == ' ' )
4896                {
4897                        status = 's';
4898                }
4899                else
4900                {
4901                        if( status == 's' ) 
4902                        {
4903                                if( *pt == '\n' || *pt == '#' ) break;
4904                                status = 'n';
4905                                t[v] = atoi( pt );
4906                                if( t[v] == 0 )
4907                                {
4908                                        fprintf( stderr, "Format error? Sequences must be specified as 1, 2, 3...\n" );
4909                                        exit( 1 );
4910                                }
4911                                if( t[v] < 0 ) *preservegaps = 1;
4912                                t[v] = abs( t[v] );
4913                                t[v] -= 1;
4914                                v++;
4915                        }
4916                }
4917                pt++;
4918        }
4919        t[v] = -1;
4920        return( v );
4921}
4922
4923static int countspace( char *s )
4924{
4925        int v = 0;
4926        char status = 's';
4927        char *pt = s;
4928        tab2space( s );
4929        while( *pt )
4930        {
4931                if( *pt == ' ' )
4932                {
4933                        status = 's';
4934                }
4935                else
4936                {
4937                        if( status == 's' ) 
4938                        {
4939                                if( *pt == '\n' || *pt == '#' ) break;
4940                                v++;
4941                                status = 'n';
4942                                if( atoi( pt ) == 0 )
4943                                {
4944                                        fprintf( stderr, "Format error? Sequences should be specified as 1, 2, 3...\n" );
4945                                        exit( 1 );
4946                                }
4947                        }
4948                }
4949                pt++;
4950        }
4951        return( v );
4952}
4953
4954
4955void readsubalignmentstable( int nseq, int **table, int *preservegaps, int *nsubpt, int *maxmempt )
4956{
4957        FILE *fp;
4958        char *line;
4959        int linelen = 1000000;
4960        int nmem;
4961        int lpos;
4962        int i, p;
4963        int *tab01;
4964
4965        line = calloc( linelen, sizeof( char ) );
4966        fp = fopen( "_subalignmentstable", "r" );
4967        if( !fp )
4968        {
4969                fprintf( stderr, "Cannot open _subalignmentstable\n" );
4970                exit( 1 );
4971        }
4972        if( table == NULL )
4973        {
4974                *nsubpt = 0;
4975                *maxmempt = 0;
4976                while( 1 )
4977                {
4978                        fgets( line, linelen-1, fp );
4979                        if( feof( fp ) ) break;
4980                        if( line[strlen(line)-1] != '\n' )
4981                        {
4982                                fprintf( stderr, "too long line? \n" );
4983                                exit( 1 );
4984                        }
4985                        if( line[0] == '#' ) continue;
4986                        if( atoi( line ) == 0 ) continue;
4987                        nmem = countspace( line );
4988                        if( nmem > *maxmempt ) *maxmempt = nmem;
4989                        (*nsubpt)++;
4990                }
4991        }
4992        else
4993        {
4994                tab01 = calloc( nseq, sizeof( int ) );
4995                for( i=0; i<nseq; i++ ) tab01[i] = 0;
4996                lpos = 0;
4997                while( 1 )
4998                {
4999                        fgets( line, linelen-1, fp );
5000                        if( feof( fp ) ) break;
5001                        if( line[strlen(line)-1] != '\n' )
5002                        {
5003                                fprintf( stderr, "too long line? \n" );
5004                                exit( 1 );
5005                        }
5006                        if( line[0] == '#' ) continue;
5007                        if( atoi( line ) == 0 ) continue;
5008                        readasubalignment( line, table[lpos], preservegaps+lpos );
5009                        for( i=0; (p=table[lpos][i])!=-1; i++ )
5010                        {
5011                                if( tab01[p] )
5012                                {
5013                                        fprintf( stderr, "\nSequence %d appears in different groups.\n", p+1 );
5014                                        fprintf( stderr, "Hierarchical grouping is not supported.\n\n" );
5015                                        exit( 1 );
5016                                }
5017                                tab01[p] = 1;
5018                                if( p > nseq-1 )
5019                                {
5020                                        fprintf( stderr, "Sequence %d does not exist in the input sequence file.\n", p+1 );
5021                                        exit( 1 );
5022                                }
5023                        }
5024                        lpos++;
5025                }
5026                free( tab01 );
5027        }
5028        fclose( fp );
5029        free( line );
5030}
5031
5032
5033void readmccaskill( FILE *fp, RNApair **pairprob, int length )
5034{
5035        char gett[1000];
5036        int *pairnum;
5037        int i;
5038        int left, right;
5039        float prob;
5040        int c;
5041
5042        pairnum = (int *)calloc( length, sizeof( int ) );
5043        for( i=0; i<length; i++ ) pairnum[i] = 0;
5044
5045        c = getc( fp );
5046        {
5047                if( c != '>' )
5048                {
5049                        fprintf( stderr, "format error in hat4 - 1\n" );
5050                        exit( 1 );
5051                }
5052        }
5053        fgets( gett, 999, fp );
5054        while( 1 )
5055        {
5056                if( feof( fp ) ) break;
5057                c = getc( fp );
5058                ungetc( c, fp );
5059                if( c == '>' || c == EOF )
5060                {
5061                        break;
5062                }
5063                fgets( gett, 999, fp );
5064//              fprintf( stderr, "gett = %s\n", gett );
5065                sscanf( gett, "%d %d %f", &left, &right, &prob );
5066
5067                if( left >= length || right >= length )
5068                {
5069                        fprintf( stderr, "format error in hat4 - 2\n" );
5070                        exit( 1 );
5071                }
5072
5073                if( prob < 0.01 ) continue; // 080607, mafft ni dake eikyou
5074
5075                if( left != right && prob > 0.0 )
5076                {
5077                        pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
5078                        pairprob[left][pairnum[left]].bestscore = prob;
5079                        pairprob[left][pairnum[left]].bestpos = right;
5080                        pairnum[left]++;
5081                        pairprob[left][pairnum[left]].bestscore = -1.0;
5082                        pairprob[left][pairnum[left]].bestpos = -1;
5083//                      fprintf( stderr, "%d-%d, %f\n", left, right, prob );
5084
5085                        pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
5086                        pairprob[right][pairnum[right]].bestscore = prob;
5087                        pairprob[right][pairnum[right]].bestpos = left;
5088                        pairnum[right]++;
5089                        pairprob[right][pairnum[right]].bestscore = -1.0;
5090                        pairprob[right][pairnum[right]].bestpos = -1;
5091//                      fprintf( stderr, "%d-%d, %f\n", right, left, prob );
5092                }
5093        }
5094        free( pairnum );
5095}
5096
5097void readpairfoldalign( FILE *fp, char *s1, char *s2, char *aln1, char *aln2, int q1, int q2, int *of1, int *of2, int sumlen )
5098{
5099        char gett[1000];
5100        int *maptoseq1;
5101        int *maptoseq2;
5102        char dumc;
5103        int dumi;
5104        char sinseq[100], sinaln[100];
5105        int posinseq, posinaln;
5106        int alnlen;
5107        int i;
5108        int pos1, pos2;
5109        char *pa1, *pa2;
5110        char qstr[1000];
5111
5112        *of1 = -1;
5113        *of2 = -1;
5114
5115        maptoseq1 = AllocateIntVec( sumlen+1 );
5116        maptoseq2 = AllocateIntVec( sumlen+1 );
5117
5118        posinaln = 0; // foldalign ga alingment wo kaesanaitok no tame.
5119
5120        while( !feof( fp ) )
5121        {
5122                fgets( gett, 999, fp );
5123                if( !strncmp( gett, "; ALIGNING", 10 ) ) break;
5124        }
5125        sprintf( qstr, "; ALIGNING            %d against %d\n", q1+1, q2+1 );
5126        if( strcmp( gett, qstr ) )
5127        {
5128                fprintf( stderr, "Error in FOLDALIGN\n" );
5129                fprintf( stderr, "qstr = %s, but gett = %s\n", qstr, gett );
5130                exit( 1 );
5131        }
5132
5133        while( !feof( fp ) )
5134        {
5135                fgets( gett, 999, fp );
5136                if( !strncmp( gett, "; --------", 10 ) ) break;
5137        }
5138
5139
5140        while( !feof( fp ) )
5141        {
5142                fgets( gett, 999, fp );
5143                if( !strncmp( gett, "; ********", 10 ) ) break;
5144//              fprintf( stderr, "gett = %s\n", gett );
5145                sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
5146                posinaln = atoi( sinaln );
5147                posinseq = atoi( sinseq );
5148//              fprintf( stderr, "posinseq = %d\n", posinseq );
5149//              fprintf( stderr, "posinaln = %d\n", posinaln );
5150                maptoseq1[posinaln-1] = posinseq-1;
5151        }
5152        alnlen = posinaln;
5153
5154        while( !feof( fp ) )
5155        {
5156                fgets( gett, 999, fp );
5157                if( !strncmp( gett, "; --------", 10 ) ) break;
5158        }
5159
5160        while( !feof( fp ) )
5161        {
5162                fgets( gett, 999, fp );
5163                if( !strncmp( gett, "; ********", 10 ) ) break;
5164//              fprintf( stderr, "gett = %s\n", gett );
5165                sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
5166                posinaln = atof( sinaln );
5167                posinseq = atof( sinseq );
5168//              fprintf( stderr, "posinseq = %d\n", posinseq );
5169//              fprintf( stderr, "posinaln = %d\n", posinaln );
5170                maptoseq2[posinaln-1] = posinseq-1;
5171        }
5172        if( alnlen != posinaln )
5173        {
5174                fprintf( stderr, "Error in foldalign?\n" );
5175                exit( 1 );
5176        }
5177
5178        pa1 = aln1;
5179        pa2 = aln2;
5180        for( i=0; i<alnlen; i++ )
5181        {
5182                pos1 = maptoseq1[i];
5183                pos2 = maptoseq2[i];
5184
5185                if( pos1 > -1 )
5186                        *pa1++ = s1[pos1];
5187                else
5188                        *pa1++ = '-';
5189
5190                if( pos2 > -1 )
5191                        *pa2++ = s2[pos2];
5192                else
5193                        *pa2++ = '-';
5194        }
5195        *pa1 = 0;
5196        *pa2 = 0;
5197
5198        *of1 = 0;
5199        for( i=0; i<alnlen; i++ )
5200        {
5201                *of1 = maptoseq1[i];
5202                if( *of1 > -1 ) break;
5203        }
5204        *of2 = 0;
5205        for( i=0; i<alnlen; i++ )
5206        {
5207                *of2 = maptoseq2[i];
5208                if( *of2 > -1 ) break;
5209        }
5210
5211//      fprintf( stderr, "*of1=%d, aln1 = :%s:\n", *of1, aln1 );
5212//      fprintf( stderr, "*of2=%d, aln2 = :%s:\n", *of2, aln2 );
5213
5214        free( maptoseq1 );
5215        free( maptoseq2 );
5216}
Note: See TracBrowser for help on using the repository browser.