binaries/src/mafft/core/io.c

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