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 == '=' || 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 == '=' || c == EOF ) && b == '\n' ) )
1189                 b = c;
1190         ungetc( c, fp );
1191 }
1192
1193 static int onlyAlpha_lower( char *str )
1194 {
1195         char tmp;
1196         char *res = str;
1197         char *bk = str;
1198
1199         while( (tmp=*str++) )
1200                 if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1201                         *res++ = tolower( tmp );
1202         *res = 0;
1203         return( res - bk );
1204 }
1205 static int onlyAlpha_upper( char *str )
1206 {
1207         char tmp;
1208         char *res = str;
1209         char *bk = str;
1210
1211         while( (tmp=*str++) )
1212                 if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1213                         *res++ = toupper( tmp );
1214         *res = 0;
1215         return( res - bk );
1216 }
1217
1218 void kake2hiku( char *str )
1219 {
1220         do
1221                 if( *str == '*' ) *str = '-';
1222         while( *str++ );
1223 }
1224
1225 char *load1SeqWithoutName_realloc( FILE *fpp )
1226 {
1227         int c, b;
1228         char *cbuf;
1229         int size = N;
1230         char *val;
1231
1232         val = malloc( (size+1) * sizeof( char ) );
1233         cbuf = val;
1234
1235         b = '\n';
1236         while( ( c = getc( fpp ) ) != EOF &&
1237           !( ( c == '>' || c == '=' || c == '(' || c == EOF ) && b == '\n' ) )
1238         {
1239                 *cbuf++ = (char)c;  /* Ä¹¤¹¤®¤Æ¤â¤·¤é¤Ê¤¤ */
1240                 if( cbuf - val == size )
1241                 {
1242                         size += N;
1243                         fprintf( stderr, "reallocating...\n" );
1244                         val = (char *)realloc( val, (size+1) * sizeof( char ) );
1245                         if( !val )
1246                         {
1247                                 fprintf( stderr, "Allocation error in load1SeqWithoutName_realloc \n" );
1248                                 exit( 1 );
1249                         }
1250                         fprintf( stderr, "done.\n" );
1251                         cbuf = val + size-N;
1252                 }
1253                 b = c;
1254         }
1255         ungetc( c, fpp );
1256         *cbuf = 0;
1257         if( dorp == 'd' )
1258                 onlyAlpha_lower( val );
1259         else
1260                 onlyAlpha_upper( val );
1261         kake2hiku( val );
1262         return( val );
1263 }
1264
1265 int load1SeqWithoutName_new( FILE *fpp, char *cbuf )
1266 {
1267         int c, b;
1268         char *bk = cbuf;
1269
1270         b = '\n';
1271         while( ( c = getc( fpp ) ) != EOF &&                    /* by T. Nishiyama */
1272           !( ( c == '>' || c == '=' || c == '(' || c == EOF ) && b == '\n' ) )
1273         {
1274                 *cbuf++ = (char)c;  /* Ä¹¤¹¤®¤Æ¤â¤·¤é¤Ê¤¤ */
1275                 b = c;
1276         }
1277         ungetc( c, fpp );
1278         *cbuf = 0;
1279         if( dorp == 'd' )
1280                 onlyAlpha_lower( bk );
1281         else
1282                 onlyAlpha_upper( bk );
1283         kake2hiku( bk );
1284         return( 0 );
1285 }
1286
1287
1288 void readDataforgaln( FILE *fp, char **name, int *nlen, char **seq )
1289 {
1290         int i;
1291         static char *tmpseq = NULL;
1292
1293 #if 0
1294         if( !tmpseq )
1295         {
1296                 tmpseq = AllocateCharVec( N );
1297         }
1298 #endif
1299
1300         rewind( fp );
1301         searchKUorWA( fp );
1302
1303         for( i=0; i<njob; i++ )
1304         {
1305                 name[i][0] = '='; getc( fp );
1306 #if 0
1307                 fgets( name[i]+1, B-2, fp );
1308                 j = strlen( name[i] );
1309                 if( name[i][j-1] != '\n' )
1310                         ErrorExit( "Too long name\n" );
1311                 name[i][j-1] = 0;
1312 #else
1313                 myfgets( name[i]+1, B-2, fp );
1314 #endif
1315 #if 0
1316                 fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1317 #endif
1318                 tmpseq = load1SeqWithoutName_realloc( fp );
1319                 strcpy( seq[i], tmpseq );
1320                 nlen[i] = strlen( seq[i] );
1321                 free( tmpseq );
1322         }
1323         if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1324 #if 0
1325         free( tmpseq );
1326 #endif
1327 }
1328
1329 void readData_varlen( FILE *fp, char **name, int *nlen, char **seq )
1330 {
1331         int i;
1332         static char *tmpseq = NULL;
1333
1334         rewind( fp );
1335         searchKUorWA( fp );
1336
1337         for( i=0; i<njob; i++ )
1338         {
1339                 name[i][0] = '='; getc( fp );
1340 #if 0
1341                 fgets( name[i]+1, B-2, fp );
1342                 j = strlen( name[i] );
1343                 if( name[i][j-1] != '\n' )
1344                         ErrorExit( "Too long name\n" );
1345                 name[i][j-1] = 0;
1346 #else
1347                 myfgets( name[i]+1, B-2, fp );
1348 #endif
1349 #if 0
1350                 fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1351 #endif
1352                 tmpseq = load1SeqWithoutName_realloc( fp );
1353                 nlen[i] = strlen( tmpseq );
1354 //              fprintf( stderr, "nlen[%d] = %d\n", i+1, nlen[i] );
1355                 seq[i] = calloc( nlen[i]+1, sizeof( char ) );
1356                 strcpy( seq[i], tmpseq );
1357                 free( tmpseq );
1358         }
1359         if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1360 #if 0
1361         free( tmpseq );
1362 #endif
1363 }
1364
1365 void readData_pointer( FILE *fp, char **name, int *nlen, char **seq )
1366 {
1367         int i;
1368         static char *tmpseq = NULL;
1369
1370 #if 0
1371         if( !tmpseq )
1372         {
1373                 tmpseq = AllocateCharVec( N );
1374         }
1375 #endif
1376
1377         rewind( fp );
1378         searchKUorWA( fp );
1379
1380         for( i=0; i<njob; i++ )
1381         {
1382                 name[i][0] = '='; getc( fp );
1383 #if 0
1384                 fgets( name[i]+1, B-2, fp );
1385                 j = strlen( name[i] );
1386                 if( name[i][j-1] != '\n' )
1387                         ErrorExit( "Too long name\n" );
1388                 name[i][j-1] = 0;
1389 #else
1390                 myfgets( name[i]+1, B-2, fp );
1391 #endif
1392 #if 0
1393                 fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1394 #endif
1395                 tmpseq = load1SeqWithoutName_realloc( fp );
1396                 strcpy( seq[i], tmpseq );
1397                 free( tmpseq );
1398                 nlen[i] = strlen( seq[i] );
1399         }
1400         if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1401 #if 0
1402         free( tmpseq );
1403 #endif
1404 }
1405
1406 void readData( FILE *fp, char name[][B], int nlen[], char **seq )
1407 {
1408         int i;
1409         static char *tmpseq = NULL;
1410
1411 #if 0
1412         if( !tmpseq )
1413         {
1414                 tmpseq = AllocateCharVec( N );
1415         }
1416 #endif
1417
1418         rewind( fp );
1419         searchKUorWA( fp );
1420
1421         for( i=0; i<njob; i++ )
1422         {
1423                 name[i][0] = '='; getc( fp );
1424 #if 0
1425                 fgets( name[i]+1, B-2, fp );
1426                 j = strlen( name[i] );
1427                 if( name[i][j-1] != '\n' )
1428                         ErrorExit( "Too long name\n" );
1429                 name[i][j-1] = 0;
1430 #else
1431                 myfgets( name[i]+1, B-2, fp );
1432 #endif
1433 #if 0
1434                 fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1435 #endif
1436                 tmpseq = load1SeqWithoutName_realloc( fp );
1437                 strcpy( seq[i], tmpseq );
1438                 nlen[i] = strlen( seq[i] );
1439                 free( tmpseq );
1440         }
1441         if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1442 #if 0
1443         free( tmpseq );
1444 #endif
1445 }
1446
1447 void cutAlignment( FILE *fp, int **regtable, char **revtable, int *outtable, char **name, char **outseq )
1448 {
1449         int i, j;
1450         int outlen;
1451         static char *tmpseq = NULL;
1452         static char *dumname = NULL;
1453         char *fs, *rs;
1454         int npos, lpos;
1455         int startpos, endpos, seqlen;
1456
1457         if( dumname == NULL )
1458         {
1459                 dumname = AllocateCharVec( N );
1460         }
1461
1462         rewind( fp );
1463         searchKUorWA( fp );
1464
1465
1466         npos = 0;
1467         for( i=0; i<njob; i++ )
1468         {
1469                 dumname[0] = '>'; getc( fp );
1470                 myfgets( dumname+1, B-1, fp );
1471                 tmpseq = load1SeqWithoutName_realloc( fp );
1472
1473                 if( outtable[i] )
1474                 {
1475 //                      putc( '>', stdout );
1476 //                      puts( dumname+1 );
1477
1478
1479                         strncat( name[npos], dumname, B-1 );
1480                         name[npos][B-1] = 0;
1481
1482                         if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1483                         seqlen = strlen( tmpseq );
1484                         lpos = 0;
1485                         for( j=0; j<5; j++ )
1486                         {
1487                                 if( regtable[0][j*2] == -1 && regtable[0][j*2+1] == -1 ) continue;
1488
1489                                 startpos = regtable[0][j*2];
1490                                 endpos   = regtable[0][j*2+1];
1491                                 if( startpos > endpos )
1492                                 {
1493                                         endpos   = regtable[0][j*2];
1494                                         startpos = regtable[0][j*2+1];
1495                                 }
1496
1497                                 if( startpos < 0 ) startpos = 0;
1498                                 if( endpos   < 0 ) endpos   = 0;
1499                                 if( endpos   >= seqlen ) endpos   = seqlen-1;
1500                                 if( startpos >= seqlen ) startpos = seqlen-1;
1501
1502 //                              fprintf( stderr, "startpos = %d, endpos = %d\n", startpos, endpos );
1503
1504                                 outlen = endpos - startpos+1;
1505                                 if( revtable[0][j] == 'f' )
1506                                 {
1507 //                                      fprintf( stderr, "regtable[%d][st] = %d\n", i, regtable[0][j*2+0] );
1508 //                                      fprintf( stderr, "regtable[%d][en] = %d\n", i, regtable[0][j*2+1] );
1509 //                                      fprintf( stderr, "outlen = %d\n", outlen );
1510 //                                      fprintf( stdout, "%.*s\n", outlen, tmpseq+regtable[0][j*2] );
1511                                         strncpy( outseq[npos] + lpos, tmpseq+startpos, outlen );
1512                                         lpos += outlen;
1513                                 }
1514                                 else
1515                                 {
1516                                         fs = AllocateCharVec( outlen+1 );
1517                                         rs = AllocateCharVec( outlen+1 );
1518
1519                                         fs[outlen] = 0;
1520                                         strncpy( fs, tmpseq+startpos, outlen );
1521                                         sreverse( rs, fs );
1522 //                                      fprintf( stdout, "%s\n", rs );
1523                                         strncpy( outseq[npos] + lpos, rs, outlen );
1524                                         lpos += outlen;
1525                                         free( fs );
1526                                         free( rs );
1527                                 }
1528                                 outseq[npos][lpos] = 0;
1529                         }
1530                         npos++;
1531                 }
1532                 free( tmpseq );
1533         }
1534 }
1535
1536 void cutData( FILE *fp, int **regtable, char **revtable, int *outtable )
1537 {
1538         int i, j;
1539         int outlen, seqlen, startpos, endpos;
1540         static char *tmpseq = NULL;
1541         static char *dumname = NULL;
1542         char *fs, *rs;
1543
1544         if( dumname == NULL )
1545         {
1546                 dumname = AllocateCharVec( N );
1547         }
1548
1549         rewind( fp );
1550         searchKUorWA( fp );
1551
1552         for( i=0; i<njob; i++ )
1553         {
1554                 dumname[0] = '='; getc( fp );
1555                 myfgets( dumname+1, B-2, fp );
1556                 tmpseq = load1SeqWithoutName_realloc( fp );
1557
1558                 if( outtable[i] )
1559                 {
1560                         gappick_samestring( tmpseq );
1561                         putc( '>', stdout );
1562                         puts( dumname+1 );
1563
1564                         seqlen = strlen( tmpseq );
1565
1566                         if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1567                         for( j=0; j<5; j++ )
1568                         {
1569                                 if( regtable[i][j*2] == -1 && regtable[i][j*2+1] == -1 ) continue;
1570
1571                                 startpos = regtable[i][j*2];
1572                                 endpos   = regtable[i][j*2+1];
1573
1574                                 if( startpos > endpos )
1575                                 {
1576                                         endpos   = regtable[i][j*2];
1577                                         startpos = regtable[i][j*2+1];
1578                                 }
1579
1580                                 if( startpos < 0 ) startpos = 0;
1581                                 if( endpos   < 0 ) endpos   = 0;
1582                                 if( endpos   >= seqlen ) endpos   = seqlen-1;
1583                                 if( startpos >= seqlen ) startpos = seqlen-1;
1584
1585                                 outlen = endpos - startpos + 1;
1586                                 if( revtable[i][j] == 'f' )
1587                                 {
1588                                         fprintf( stderr, "startpos = %d\n", startpos );
1589                                         fprintf( stderr, "endpos   = %d\n", endpos );
1590                                         fprintf( stderr, "outlen = %d\n", outlen );
1591                                         fprintf( stdout, "%.*s\n", outlen, tmpseq+startpos );
1592                                 }
1593                                 else
1594                                 {
1595                                         fs = AllocateCharVec( outlen+1 );
1596                                         rs = AllocateCharVec( outlen+1 );
1597
1598                                         fs[outlen] = 0;
1599                                         strncpy( fs, tmpseq+startpos, outlen );
1600                                         sreverse( rs, fs );
1601                                         fprintf( stdout, "%s\n", rs );
1602                                         free( fs );
1603                                         free( rs );
1604                                 }
1605                         }
1606                 }
1607                 free( tmpseq );
1608         }
1609 }
1610
1611 void catData( FILE *fp )
1612 {
1613         int i;
1614         static char *tmpseq = NULL;
1615         static char *dumname = NULL;
1616
1617         if( dumname == NULL )
1618         {
1619                 dumname = AllocateCharVec( N );
1620         }
1621
1622         rewind( fp );
1623         searchKUorWA( fp );
1624
1625         for( i=0; i<njob; i++ )
1626         {
1627                 dumname[0] = '='; getc( fp );
1628                 myfgets( dumname+1, B-2, fp );
1629                 putc( '>', stdout );
1630                 puts( dumname+1 );
1631                 tmpseq = load1SeqWithoutName_realloc( fp );
1632                 if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1633                 puts( tmpseq );
1634                 free( tmpseq );
1635         }
1636 }
1637
1638 int countATGC( char *s, int *total )
1639 {
1640         int nATGC;
1641         int nChar;
1642         char c;
1643         nATGC = nChar = 0;
1644
1645         if( *s == 0 )
1646         {
1647                 total = 0;
1648                 return( 0 );
1649         }
1650
1651         do
1652         {
1653                 c = tolower( *s );
1654                 if( isalpha( c ) )
1655                 {
1656                         nChar++;
1657                         if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1658                                 nATGC++;
1659                 }
1660         }
1661         while( *++s );
1662
1663         *total = nChar;
1664         return( nATGC );
1665 }
1666
1667 double countATGCbk( char *s )
1668 {
1669         int nATGC;
1670         int nChar;
1671         char c;
1672         nATGC = nChar = 0;
1673
1674         do
1675         {
1676                 c = tolower( *s );
1677                 if( isalpha( c ) )
1678                 {
1679                         nChar++;
1680                         if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1681                                 nATGC++;
1682                 }
1683         }
1684         while( *++s );
1685         return( (double)nATGC / nChar );
1686 }
1687
1688
1689 int countalpha( char *seq )
1690 {
1691         int val = 0;
1692         while( *seq )
1693                 if( isalpha( *seq++ ) ) val++;
1694         return( val );
1695 }
1696
1697 void getnumlen_nogap( FILE *fp, int *nlenminpt )
1698 {
1699         int total;
1700         int nsite;
1701         int atgcnum;
1702         int i, tmp;
1703         char *tmpseq, *tmpname;
1704         double atgcfreq;
1705         tmpname = AllocateCharVec( N );
1706         njob = countKUorWA( fp );
1707         searchKUorWA( fp );
1708         nlenmax = 0;
1709         *nlenminpt = 99999999;
1710         atgcnum = 0;
1711         total = 0;
1712         for( i=0; i<njob; i++ )
1713         {
1714                 myfgets( tmpname, N-1, fp );
1715                 tmpseq = load1SeqWithoutName_realloc( fp );
1716                 tmp = countalpha( tmpseq );
1717                 if( tmp > nlenmax ) nlenmax  = tmp;
1718                 if( tmp < *nlenminpt ) *nlenminpt  = tmp;
1719                 atgcnum += countATGC( tmpseq, &nsite );
1720                 total += nsite;
1721                 free( tmpseq );
1722         }
1723         free( tmpname );
1724         atgcfreq = (double)atgcnum / total;
1725         fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1726         if( dorp == NOTSPECIFIED )
1727         {
1728                 if( atgcfreq > 0.75 )
1729                 {
1730                         dorp = 'd';
1731                         upperCase = -1;
1732                 }
1733                 else
1734                 {
1735                         dorp = 'p';
1736                         upperCase = 0;
1737                 }
1738         }
1739 }
1740
1741
1742 void getnumlen_nogap_outallreg( FILE *fp, int *nlenminpt )
1743 {
1744         int total;
1745         int nsite;
1746         int atgcnum;
1747         int i, tmp;
1748         char *tmpseq, *tmpname;
1749         double atgcfreq;
1750         tmpname = AllocateCharVec( N );
1751         njob = countKUorWA( fp );
1752         searchKUorWA( fp );
1753         nlenmax = 0;
1754         *nlenminpt = 99999999;
1755         atgcnum = 0;
1756         total = 0;
1757         for( i=0; i<njob; i++ )
1758         {
1759                 myfgets( tmpname, N-1, fp );
1760                 fprintf( stdout, "%s\n", tmpname );
1761                 tmpseq = load1SeqWithoutName_realloc( fp );
1762                 tmp = countalpha( tmpseq );
1763                 fprintf( stdout, "%d\n", tmp );
1764                 if( tmp > nlenmax ) nlenmax  = tmp;
1765                 if( tmp < *nlenminpt ) *nlenminpt  = tmp;
1766                 atgcnum += countATGC( tmpseq, &nsite );
1767                 total += nsite;
1768                 free( tmpseq );
1769         }
1770         free( tmpname );
1771         atgcfreq = (double)atgcnum / total;
1772         fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1773         if( dorp == NOTSPECIFIED )
1774         {
1775                 if( atgcfreq > 0.75 )
1776                 {
1777                         dorp = 'd';
1778                         upperCase = -1;
1779                 }
1780                 else
1781                 {
1782                         dorp = 'p';
1783                         upperCase = 0;
1784                 }
1785         }
1786 }
1787
1788 void getnumlen_nogap_outallreg_web( FILE *fp, FILE *ofp, int *nlenminpt, int *isalignedpt )
1789 {
1790         int total;
1791         int nsite;
1792         int atgcnum;
1793         int alnlen, alnlen_prev;
1794         int i, tmp;
1795         char *tmpseq, *tmpname;
1796         double atgcfreq;
1797         tmpname = AllocateCharVec( N );
1798         njob = countKUorWA( fp );
1799         searchKUorWA( fp );
1800         nlenmax = 0;
1801         *nlenminpt = 99999999;
1802         atgcnum = 0;
1803         total = 0;
1804         alnlen_prev = -1;
1805         *isalignedpt = 1;
1806         for( i=0; i<njob; i++ )
1807         {
1808                 myfgets( tmpname, N-1, fp );
1809 //              fprintf( stdout, "%s\n", tmpname );
1810                 tmpseq = load1SeqWithoutName_realloc( fp );
1811                 tmp = countalpha( tmpseq );
1812 //              fprintf( stdout, "%d\n", tmp );
1813                 if( tmp > nlenmax ) nlenmax  = tmp;
1814                 if( tmp < *nlenminpt ) *nlenminpt  = tmp;
1815                 atgcnum += countATGC( tmpseq, &nsite );
1816                 total += nsite;
1817
1818                 alnlen = strlen( tmpseq );
1819 //              fprintf( stdout, "##### alnlen, alnlen_prev = %d, %d\n", alnlen, alnlen_prev );
1820                 if( i>0 && alnlen_prev != alnlen ) *isalignedpt = 0;
1821                 alnlen_prev = alnlen;
1822
1823                 free( tmpseq );
1824                 atgcfreq = (double)atgcnum / total;
1825 //              fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1826 //              if( dorp == NOTSPECIFIED ) // you kentou
1827                 {
1828                         if( atgcfreq > 0.75 )
1829                         {
1830                                 dorp = 'd';
1831                                 upperCase = -1;
1832                         }
1833                         else
1834                         {
1835                                 dorp = 'p';
1836                                 upperCase = 0;
1837                         }
1838                 }
1839
1840
1841
1842
1843                 fprintf( ofp, " <label for='s%d'><input type='checkbox' id='s%d' name='s%d' checked>%s</label>\n", i, i, i, tmpname );
1844                 fprintf( ofp, "<span id='t%d-0' style='display:none'>", i );
1845                 fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"t%d\")'>+reg</a>", i );
1846                 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 );
1847                 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 );
1848                 fprintf( ofp, "  Sequence Length:<input type='text' name='l%d' size='8' value='%d' readonly='readonly'>", i, tmp );
1849                 fprintf( ofp, "\n</span>" );
1850                 fprintf( ofp, "<span id='t%d-1' style='display:none'>", i );
1851                 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 );
1852                 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 );
1853                 fprintf( ofp, "\n</span>" );
1854                 fprintf( ofp, "<span id='t%d-2' style='display:none'>", i );
1855                 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 );
1856                 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 );
1857                 fprintf( ofp, "\n</span>" );
1858                 fprintf( ofp, "<span id='t%d-3' style='display:none'>", i );
1859                 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 );
1860                 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 );
1861                 fprintf( ofp, "\n</span>" );
1862                 fprintf( ofp, "<span id='t%d-4' style='display:none'>", i );
1863                 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 );
1864                 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 );
1865                 fprintf( ofp, "\n</span>" );
1866         }
1867         free( tmpname );
1868         atgcfreq = (double)atgcnum / total;
1869         fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1870 //      if( dorp == NOTSPECIFIED ) // you kentou
1871         {
1872                 if( atgcfreq > 0.75 )
1873                 {
1874                         dorp = 'd';
1875                         upperCase = -1;
1876                 }
1877                 else
1878                 {
1879                         dorp = 'p';
1880                         upperCase = 0;
1881                 }
1882         }
1883         if( *isalignedpt )
1884         {
1885                 fprintf( ofp, "\n" );
1886                 fprintf( ofp, "<span id='tall-0' style='display:none'>" );
1887                 fprintf( ofp, "Cut the alignment\n" );
1888                 fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"tall\")'>+reg</a>" );
1889                 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 );
1890                 if( dorp == 'd' ) fprintf( ofp, " <label for='rall-0'><input type='checkbox' name='rall-0' id='rall-0'>Reverse</label>" );
1891                 fprintf( ofp, "  Alignment length:<input type='text' name='lall' size='8' value='%d' readonly='readonly'>", alnlen );
1892                 fprintf( ofp, "\n</span>" );
1893                 fprintf( ofp, "<span id='tall-1' style='display:none'>" );
1894                 fprintf( ofp, "      Begin:<input type='text' name='ball-1' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
1895                 if( dorp == 'd' ) fprintf( ofp, " <label for='rall-1'><input type='checkbox' name='rall-1' id='rall-1'>Reverse</label>" );
1896                 fprintf( ofp, "\n</span>" );
1897                 fprintf( ofp, "<span id='tall-2' style='display:none'>" );
1898                 fprintf( ofp, "      Begin:<input type='text' name='ball-2' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
1899                 if( dorp == 'd' ) fprintf( ofp, " <label for='rall-2'><input type='checkbox' name='rall-2' id='rall-2'>Reverse</label>" );
1900                 fprintf( ofp, "\n</span>" );
1901                 fprintf( ofp, "<span id='tall-3' style='display:none'>" );
1902                 fprintf( ofp, "      Begin:<input type='text' name='ball-3' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
1903                 if( dorp == 'd' ) fprintf( ofp, " <label for='rall-3'><input type='checkbox' name='rall-3' id='rall-3'>Reverse</label>" );
1904                 fprintf( ofp, "\n</span>" );
1905                 fprintf( ofp, "<span id='tall-4' style='display:none'>" );
1906                 fprintf( ofp, "      Begin:<input type='text' name='ball-4' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
1907                 if( dorp == 'd' ) fprintf( ofp, " <label for='rall-4'><input type='checkbox' name='rall-4' id='rall-4'>Reverse</label>" );
1908                 fprintf( ofp, "\n</span>" );
1909         }
1910
1911 }
1912
1913 void getnumlen( FILE *fp )
1914 {
1915         int total;
1916         int nsite;
1917         int atgcnum;
1918         int i, tmp;
1919         char *tmpseq;
1920         char *tmpname;
1921         double atgcfreq;
1922         tmpname = AllocateCharVec( N );
1923         njob = countKUorWA( fp );
1924         searchKUorWA( fp );
1925         nlenmax = 0;
1926         atgcnum = 0;
1927         total = 0;
1928         for( i=0; i<njob; i++ )
1929         {
1930                 myfgets( tmpname, N-1, fp );
1931                 tmpseq = load1SeqWithoutName_realloc( fp );
1932                 tmp = strlen( tmpseq );
1933                 if( tmp > nlenmax ) nlenmax  = tmp;
1934                 atgcnum += countATGC( tmpseq, &nsite );
1935                 total += nsite;
1936 //              fprintf( stderr, "##### total = %d\n", total );
1937                 free( tmpseq );
1938         }
1939
1940         atgcfreq = (double)atgcnum / total;
1941 //      fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1942         if( dorp == NOTSPECIFIED )
1943         {
1944                 if( atgcfreq > 0.75 )
1945                 {
1946                         dorp = 'd';
1947                         upperCase = -1;
1948                 }
1949                 else
1950                 {
1951                         dorp = 'p';
1952                         upperCase = 0;
1953                 }
1954         }
1955         free( tmpname );
1956 }
1957
1958
1959
1960 void WriteGapFill( FILE *fp, int locnjob, char name[][B], int nlen[M], char **aseq )
1961 {
1962         static char b[N];
1963         int i, j;
1964         int nalen[M];
1965         static char gap[N];
1966         static char buff[N];
1967
1968 #if IODEBUG
1969         fprintf( stderr, "IMAKARA KAKU\n" );
1970 #endif
1971         nlenmax = 0;
1972         for( i=0; i<locnjob; i++ )
1973         {
1974                 int len = strlen( aseq[i] );
1975                 if( nlenmax < len ) nlenmax = len;
1976         }
1977
1978         for( i=0; i<nlenmax; i++ ) gap[i] = '-';
1979         gap[nlenmax] = 0;
1980
1981         fprintf( fp, "%5d", locnjob );
1982         fprintf( fp, "\n" );
1983
1984         for( i=0; i<locnjob; i++ )
1985         {
1986                 strcpy( buff, aseq[i] );
1987                 strncat( buff, gap, nlenmax-strlen( aseq[i] ) );
1988                 buff[nlenmax] = 0;
1989                 nalen[i] = strlen( buff );
1990                 fprintf( fp, "%s\n", name[i] );
1991                 fprintf( fp, "%5d\n", nalen[i] );
1992                 for( j=0; j<nalen[i]; j=j+C )
1993                 {
1994                         strncpy_caseC( b, buff+j, C ); b[C] = 0;
1995                         fprintf( fp, "%s\n",b );
1996                 }
1997         }
1998 #if DEBUG
1999         fprintf( stderr, "nalen[0] = %d\n", nalen[0] );
2000 #endif
2001 #if IODEBUG
2002         fprintf( stderr, "KAKIOWATTA\n" );
2003 #endif
2004 }
2005
2006 void writeDataforgaln( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2007 {
2008         int i, j;
2009         int nalen;
2010
2011         for( i=0; i<locnjob; i++ )
2012         {
2013                 nalen = strlen( aseq[i] );
2014                 fprintf( fp, ">%s\n", name[i]+1 );
2015                 for( j=0; j<nalen; j=j+C )
2016                 {
2017 #if 0
2018                         strncpy( b, aseq[i]+j, C ); b[C] = 0;
2019                         fprintf( fp, "%s\n",b );
2020 #else
2021                         fprintf( fp, "%.*s\n", C, aseq[i]+j );
2022 #endif
2023                 }
2024         }
2025 }
2026
2027 void writeData_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2028 {
2029         int i, j;
2030         int nalen;
2031
2032         for( i=0; i<locnjob; i++ )
2033         {
2034 #if DEBUG
2035                 fprintf( stderr, "i = %d in writeData\n", i );
2036 #endif
2037                 nalen = strlen( aseq[i] );
2038                 fprintf( fp, ">%s\n", name[i]+1 );
2039                 for( j=0; j<nalen; j=j+C )
2040                 {
2041 #if 0
2042                         strncpy( b, aseq[i]+j, C ); b[C] = 0;
2043                         fprintf( fp, "%s\n",b );
2044 #else
2045                         fprintf( fp, "%.*s\n", C, aseq[i]+j );
2046 #endif
2047                 }
2048         }
2049 }
2050
2051 void writeData( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq )
2052 {
2053         int i, j;
2054         int nalen;
2055
2056         for( i=0; i<locnjob; i++ )
2057         {
2058 #if DEBUG
2059                 fprintf( stderr, "i = %d in writeData\n", i );
2060 #endif
2061                 nalen = strlen( aseq[i] );
2062                 fprintf( fp, ">%s\n", name[i]+1 );
2063                 for( j=0; j<nalen; j=j+C )
2064                 {
2065 #if 0
2066                         strncpy( b, aseq[i]+j, C ); b[C] = 0;
2067                         fprintf( fp, "%s\n",b );
2068 #else
2069                         fprintf( fp, "%.*s\n", C, aseq[i]+j );
2070 #endif
2071                 }
2072         }
2073 }
2074
2075
2076 void write1seq( FILE *fp, char *aseq )
2077 {
2078         int j;
2079         int nalen;
2080
2081         nalen = strlen( aseq );
2082         for( j=0; j<nalen; j=j+C )
2083                 fprintf( fp, "%.*s\n", C, aseq+j );
2084 }
2085
2086
2087
2088 void readhat2_floathalf_pointer( FILE *fp, int nseq, char **name, float **mtx )
2089 {
2090     int i, j, nseq0;
2091     char b[B];
2092
2093     fgets( b, B, fp );
2094     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2095     fgets( b, B, fp );
2096     for( i=0; i<nseq; i++ )
2097     {
2098 #if 0
2099         getaline_fp_eof( b, B, fp );
2100 #else
2101                 myfgets( b, B-2, fp );
2102 #endif
2103 #if 0
2104                 j = MIN( strlen( b+6 ), 10 );
2105         if( strncmp( name[i], b+6 , j ) )
2106                 {
2107                         fprintf( stderr, "Error in hat2\n" );
2108                         fprintf( stderr, "%s != %s\n", b, name[i] );
2109                         exit( 1 );
2110                 }
2111 #endif
2112     }
2113     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2114     {
2115         mtx[i][j-i] = ( input_new( fp, D ) );
2116     }
2117 }
2118 void readhat2_floathalf( FILE *fp, int nseq, char name[M][B], float **mtx )
2119 {
2120     int i, j, nseq0;
2121     char b[B];
2122
2123     fgets( b, B, fp );
2124     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2125     fgets( b, B, fp );
2126     for( i=0; i<nseq; i++ )
2127     {
2128 #if 0
2129         getaline_fp_eof( b, B, fp );
2130 #else
2131                 myfgets( b, B-2, fp );
2132 #endif
2133 #if 0
2134                 j = MIN( strlen( b+6 ), 10 );
2135         if( strncmp( name[i], b+6 , j ) )
2136                 {
2137                         fprintf( stderr, "Error in hat2\n" );
2138                         fprintf( stderr, "%s != %s\n", b, name[i] );
2139                         exit( 1 );
2140                 }
2141 #endif
2142     }
2143     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2144     {
2145         mtx[i][j-i] = ( input_new( fp, D ) );
2146     }
2147 }
2148 void readhat2_float( FILE *fp, int nseq, char name[M][B], float **mtx )
2149 {
2150     int i, j, nseq0;
2151     char b[B];
2152
2153     fgets( b, B, fp );
2154     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2155     fgets( b, B, fp );
2156     for( i=0; i<nseq; i++ )
2157     {
2158 #if 0
2159         getaline_fp_eof( b, B, fp );
2160 #else
2161                 myfgets( b, B-2, fp );
2162 #endif
2163 #if 0
2164                 j = MIN( strlen( b+6 ), 10 );
2165         if( strncmp( name[i], b+6 , j ) )
2166                 {
2167                         fprintf( stderr, "Error in hat2\n" );
2168                         fprintf( stderr, "%s != %s\n", b, name[i] );
2169                         exit( 1 );
2170                 }
2171 #endif
2172     }
2173     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2174     {
2175         mtx[i][j] = ( input_new( fp, D ) );
2176     }
2177 }
2178 void readhat2_int( FILE *fp, int nseq, char name[M][B], int **mtx )
2179 {
2180     int i, j, nseq0;
2181     char b[B];
2182
2183     fgets( b, B, fp );
2184     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2185     fgets( b, B, fp );
2186     for( i=0; i<nseq; i++ )
2187     {
2188 #if 0
2189         getaline_fp_eof( b, B, fp );
2190 #else
2191                 myfgets( b, B-2, fp );
2192 #endif
2193 #if 0
2194                 j = MIN( strlen( b+6 ), 10 );
2195         if( strncmp( name[i], b+6 , j ) )
2196                 {
2197                         fprintf( stderr, "Error in hat2\n" );
2198                         fprintf( stderr, "%s != %s\n", b, name[i] );
2199                         exit( 1 );
2200                 }
2201 #endif
2202     }
2203     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2204     {
2205         mtx[i][j] = (int)( input_new( fp, D ) * INTMTXSCALE + 0.5 );
2206     }
2207 }
2208 void readhat2( FILE *fp, int nseq, char name[M][B], double **mtx )
2209 {
2210     int i, j, nseq0;
2211     char b[B];
2212
2213     fgets( b, B, fp );
2214     fgets( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2215     fgets( b, B, fp );
2216     for( i=0; i<nseq; i++ )
2217     {
2218 #if 0
2219         getaline_fp_eof( b, B, fp );
2220 #else
2221                 myfgets( b, B-2, fp );
2222 #endif
2223 #if 0
2224                 j = MIN( strlen( b+6 ), 10 );
2225         if( strncmp( name[i], b+6 , j ) )
2226                 {
2227                         fprintf( stderr, "Error in hat2\n" );
2228                         fprintf( stderr, "%s != %s\n", b, name[i] );
2229                         exit( 1 );
2230                 }
2231 #endif
2232     }
2233     for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2234     {
2235         mtx[i][j] = (double)input_new( fp, D);
2236     }
2237 }
2238
2239 void WriteFloatHat2_pointer_halfmtx( FILE *hat2p, int locnjob, char **name, float **mtx )
2240 {
2241         int i, j, ijsa;
2242         double max = 0.0;
2243         for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2244
2245         fprintf( hat2p, "%5d\n", 1 );
2246         fprintf( hat2p, "%5d\n", locnjob );
2247         fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2248
2249         for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2250         for( i=0; i<locnjob; i++ )
2251         {
2252                 for( j=i+1; j<njob; j++ )
2253                 {
2254                         fprintf( hat2p, "%#6.3f", mtx[i][j-i] );
2255                         ijsa = j-i;
2256                         if( ijsa % 12 == 0 || ijsa == locnjob-i-1 ) fprintf( hat2p, "\n" );
2257                 }
2258         }
2259 }
2260
2261 void WriteFloatHat2_pointer( FILE *hat2p, int locnjob, char **name, float **mtx )
2262 {
2263         int i, j;
2264         double max = 0.0;
2265         for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2266
2267         fprintf( hat2p, "%5d\n", 1 );
2268         fprintf( hat2p, "%5d\n", locnjob );
2269         fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2270
2271         for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2272         for( i=0; i<locnjob; i++ )
2273         {
2274                 for( j=1; j<locnjob-i; j++ )
2275                 {
2276                         fprintf( hat2p, "%#6.3f", mtx[i][j] );
2277                         if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2278                 }
2279         }
2280 }
2281
2282 void WriteFloatHat2( FILE *hat2p, int locnjob, char name[M][B], float **mtx )
2283 {
2284         int i, j;
2285         double max = 0.0;
2286         for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2287
2288         fprintf( hat2p, "%5d\n", 1 );
2289         fprintf( hat2p, "%5d\n", locnjob );
2290         fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2291
2292         for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2293         for( i=0; i<locnjob; i++ )
2294         {
2295                 for( j=1; j<locnjob-i; j++ )
2296                 {
2297                         fprintf( hat2p, "%#6.3f", mtx[i][j] );
2298                         if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2299                 }
2300         }
2301 }
2302
2303 void WriteHat2_int( FILE *hat2p, int locnjob, char name[M][B], int **mtx )
2304 {
2305         int i, j;
2306         double max = 0.0;
2307         for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2308         max /= INTMTXSCALE;
2309
2310         fprintf( hat2p, "%5d\n", 1 );
2311         fprintf( hat2p, "%5d\n", locnjob );
2312         fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2313
2314         for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2315         for( i=0; i<locnjob-1; i++ )
2316         {
2317                 for( j=i+1; j<locnjob; j++ )
2318                 {
2319                         fprintf( hat2p, "%#6.3f", (float)mtx[i][j] / INTMTXSCALE );
2320                         if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2321                 }
2322         }
2323 }
2324 void WriteHat2( FILE *hat2p, int locnjob, char name[M][B], double **mtx )
2325 {
2326         int i, j;
2327         double max = 0.0;
2328         for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2329
2330         fprintf( hat2p, "%5d\n", 1 );
2331         fprintf( hat2p, "%5d\n", locnjob );
2332         fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2333
2334         for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2335         for( i=0; i<locnjob-1; i++ )
2336         {
2337                 for( j=i+1; j<locnjob; j++ )
2338                 {
2339                         fprintf( hat2p, "%#6.3f", mtx[i][j] );
2340                         if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2341                 }
2342         }
2343 }
2344
2345 #if 0
2346 void WriteHat2plain( FILE *hat2p, int locnjob, double **mtx )
2347 {
2348         int i, j, ilim;
2349
2350         ilim = locnjob-1;
2351         for( i=0; i<ilim; i++ )
2352         {
2353                 fprintf( hat2p, "%d-%d d=%.3f\n", i+1, i+1, 0.0 );
2354                 for( j=i+1; j<locnjob; j++ )
2355                 {
2356                         fprintf( hat2p, "%d-%d d=%.3f\n", i+1, j+1, mtx[i][j] );
2357                 }
2358         }
2359 }
2360 #endif
2361
2362 int ReadFasta_sub( FILE *fp, double *dis, int nseq, char name[M][B] )
2363 {
2364     int i, count=0;
2365     char b[B];
2366     int junban[M];
2367
2368     count = 0;
2369     for( i=0; i<10000000 && count<nseq; i++ )
2370     {
2371         fgets( b, B-1, fp );
2372         if( !strncmp( "+==========+", b, 12 ) )
2373         {
2374             junban[count] = atoi( b+12 );
2375             count++;
2376         }
2377     }
2378
2379         for( i=0; i<nseq; i++ ) dis[i] = 0.0;
2380     count = 0;
2381     for( i=0; i<100000 && count<nseq; i++ )
2382     {
2383                 if( fgets( b, B-1, fp ) ) break;
2384         if( !strncmp( name[junban[count]], b, 20  ) )
2385         {
2386             fgets( b, B-1, fp );
2387             dis[junban[count]] = atof( b );
2388             count++;
2389         }
2390     }
2391     return 0;
2392 }
2393
2394
2395 int ReadSsearch( FILE *fp, double *dis, int nseq, char name[M][B] )
2396 {
2397     int i, count=0;
2398     char b[B];
2399     int junban[M];
2400         int opt;
2401
2402     count = 0;
2403     for( i=0; i<10000000 && count<nseq; i++ )
2404     {
2405         fgets( b, B-1, fp );
2406         if( !strncmp( "+==========+", b, 12 ) )
2407         {
2408             junban[count] = atoi( b+12 );
2409                         sscanf( b+75, "%d", &opt );
2410             dis[junban[count]] = (double)opt;
2411             count++;
2412         }
2413     }
2414
2415 /*
2416     count = 0;
2417     for( i=0; i<100000 && count<nseq; i++ )
2418     {
2419         fgets( b, B-1, fp );
2420         if( !strncmp( name[junban[count]], b, 20  ) )
2421         {
2422             dis[junban[count]] = atof( b+65 );
2423             count++;
2424         }
2425     }
2426 */
2427     return 0;
2428 }
2429
2430 int ReadBlastm7_avscore( FILE *fp, double *dis, int nin )
2431 {
2432     int count=0;
2433     char b[B];
2434         char *pt;
2435     int *junban;
2436         double score, sumscore;
2437         double len, sumlen;
2438         int qstart, qend, tstart, tend;
2439         double scorepersite;
2440         static char qal[N], tal[N], al[N];
2441         int nlocalhom;
2442
2443         junban = calloc( nin, sizeof( int ) );
2444
2445         count = 0;
2446         sumscore = 0.0;
2447         sumlen = 0.0;
2448         score = 0.0;
2449         len = 0.0;
2450         scorepersite = 0.0; // by D.Mathog, a guess
2451     while( 1 )
2452         {
2453
2454                 if( feof( fp ) ) break;
2455
2456                 while( fgets( b, B-1, fp ) )
2457                 {
2458                         if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
2459                 }
2460
2461                 if( !strncmp( "          <Hit_def>", b, 19 ) )
2462                 {
2463                         junban[count] = atoi( b+31 );
2464                         nlocalhom = 0;
2465                 }
2466
2467
2468                 while( fgets( b, B-1, fp ) )
2469                         if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
2470                 pt = b + 25;
2471                 score = atof( pt );
2472                 sumscore += score;
2473
2474
2475                 while( fgets( b, B-1, fp ) )
2476                         if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
2477                 pt = b + 30;
2478                 qstart = atoi( pt ) - 1;
2479
2480
2481                 while( fgets( b, B-1, fp ) )
2482                         if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
2483                 pt = b + 28;
2484                 qend = atoi( pt ) - 1;
2485
2486
2487                 while( fgets( b, B-1, fp ) )
2488                         if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
2489                 pt = b + 28;
2490                 tstart = atoi( pt ) - 1;
2491
2492
2493                 while( fgets( b, B-1, fp ) )
2494                         if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
2495                 pt = b + 26;
2496                 tend = atoi( pt ) - 1;
2497
2498
2499                 while( fgets( b, B-1, fp ) )
2500                         if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
2501                 pt = b + 29;
2502                 len = atoi( pt );
2503                 sumlen += len;
2504
2505
2506                 while( fgets( al, N-100, fp ) )
2507                         if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
2508
2509                 strcpy( qal, al+24 );
2510                 pt = qal;
2511                 while( *++pt != '<' )
2512                         ;
2513                 *pt = 0;
2514
2515
2516                 while( fgets( al, N-100, fp ) )
2517                         if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
2518
2519                 strcpy( tal, al+24 );
2520                 pt = tal;
2521                 while( *++pt != '<' )
2522                         ;
2523                 *pt = 0;
2524
2525
2526 //              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 );
2527
2528
2529                 while( fgets( b, B-1, fp ) )
2530                         if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
2531
2532
2533                 fgets( b, B-1, fp );
2534
2535
2536                 if( !strncmp( "          </Hit_hsps>", b, 21 ) )
2537                 {
2538                         dis[junban[count++]] = sumscore;
2539                         sumscore = 0.0;
2540                         fgets( b, B-1, fp );
2541                         fgets( b, B-1, fp );
2542                         scorepersite = sumscore / sumlen;
2543                         if( scorepersite != (int)scorepersite )
2544                         {
2545                                 fprintf( stderr, "ERROR! sumscore=%f, sumlen=%f, and scorepersite=%f\n", sumscore, sumlen, scorepersite );
2546                                 exit( 1 );
2547                         }
2548
2549                         if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
2550                 }
2551         }
2552
2553         free( junban );
2554
2555     return (int)scorepersite;
2556 }
2557 int ReadBlastm7_scoreonly( FILE *fp, double *dis, int nin )
2558 {
2559     int count=0;
2560     char b[B];
2561         char *pt;
2562     int *junban;
2563         int overlapaa;
2564         double score, sumscore;
2565         int qstart, qend, tstart, tend;
2566         static char qal[N], tal[N], al[N];
2567         int nlocalhom;
2568
2569         junban = calloc( nin, sizeof( int ) );
2570
2571         count = 0;
2572         sumscore = 0.0;
2573         score = 0.0;
2574     while( 1 )
2575         {
2576
2577                 if( feof( fp ) ) break;
2578
2579                 while( fgets( b, B-1, fp ) )
2580                 {
2581                         if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
2582                 }
2583
2584                 if( !strncmp( "          <Hit_def>", b, 19 ) )
2585                 {
2586                         junban[count] = atoi( b+31 );
2587                         nlocalhom = 0;
2588                 }
2589
2590
2591                 while( fgets( b, B-1, fp ) )
2592                         if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
2593                 pt = b + 25;
2594                 score = atof( pt );
2595                 sumscore += score;
2596
2597
2598                 while( fgets( b, B-1, fp ) )
2599                         if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
2600                 pt = b + 30;
2601                 qstart = atoi( pt ) - 1;
2602
2603
2604                 while( fgets( b, B-1, fp ) )
2605                         if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
2606                 pt = b + 28;
2607                 qend = atoi( pt ) - 1;
2608
2609
2610                 while( fgets( b, B-1, fp ) )
2611                         if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
2612                 pt = b + 28;
2613                 tstart = atoi( pt ) - 1;
2614
2615
2616                 while( fgets( b, B-1, fp ) )
2617                         if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
2618                 pt = b + 26;
2619                 tend = atoi( pt ) - 1;
2620
2621
2622                 while( fgets( b, B-1, fp ) )
2623                         if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
2624                 pt = b + 29;
2625                 overlapaa = atoi( pt );
2626
2627
2628                 while( fgets( al, N-100, fp ) )
2629                         if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
2630
2631                 strcpy( qal, al+24 );
2632                 pt = qal;
2633                 while( *++pt != '<' )
2634                         ;
2635                 *pt = 0;
2636
2637
2638                 while( fgets( al, N-100, fp ) )
2639                         if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
2640
2641                 strcpy( tal, al+24 );
2642                 pt = tal;
2643                 while( *++pt != '<' )
2644                         ;
2645                 *pt = 0;
2646
2647
2648 //              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 );
2649
2650 //              nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
2651
2652                 while( fgets( b, B-1, fp ) )
2653                         if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
2654
2655
2656                 fgets( b, B-1, fp );
2657
2658
2659                 if( !strncmp( "          </Hit_hsps>", b, 21 ) )
2660                 {
2661                         dis[junban[count++]] = sumscore;
2662                         sumscore = 0.0;
2663                         fgets( b, B-1, fp );
2664                         fgets( b, B-1, fp );
2665                         if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
2666                 }
2667         }
2668
2669         free( junban );
2670
2671     return count;
2672 }
2673
2674 int ReadBlastm7( FILE *fp, double *dis, int qmem, char name[M][B], LocalHom *localhomlist )
2675 {
2676     int count=0;
2677     char b[B];
2678         char *pt;
2679     static int junban[M];
2680         int overlapaa;
2681         double score, sumscore;
2682         int qstart, qend, tstart, tend;
2683         static char qal[N], tal[N], al[N];
2684         int nlocalhom;
2685
2686
2687
2688         count = 0;
2689         sumscore = 0.0;
2690         score = 0.0;
2691         nlocalhom = 0;
2692     while( 1 )
2693         {
2694
2695                 if( feof( fp ) ) break;
2696
2697                 while( fgets( b, B-1, fp ) )
2698                 {
2699                         if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
2700                 }
2701
2702                 if( !strncmp( "          <Hit_def>", b, 19 ) )
2703                 {
2704                         junban[count] = atoi( b+31 );
2705                         nlocalhom = 0;
2706                 }
2707
2708
2709                 while( fgets( b, B-1, fp ) )
2710                         if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
2711                 pt = b + 25;
2712                 score = atof( pt );
2713                 sumscore += score;
2714
2715
2716                 while( fgets( b, B-1, fp ) )
2717                         if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
2718                 pt = b + 30;
2719                 qstart = atoi( pt ) - 1;
2720
2721
2722                 while( fgets( b, B-1, fp ) )
2723                         if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
2724                 pt = b + 28;
2725                 qend = atoi( pt ) - 1;
2726
2727
2728                 while( fgets( b, B-1, fp ) )
2729                         if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
2730                 pt = b + 28;
2731                 tstart = atoi( pt ) - 1;
2732
2733
2734                 while( fgets( b, B-1, fp ) )
2735                         if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
2736                 pt = b + 26;
2737                 tend = atoi( pt ) - 1;
2738
2739
2740                 while( fgets( b, B-1, fp ) )
2741                         if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
2742                 pt = b + 29;
2743                 overlapaa = atoi( pt );
2744
2745
2746                 while( fgets( al, N-100, fp ) )
2747                         if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
2748
2749                 strcpy( qal, al+24 );
2750                 pt = qal;
2751                 while( *++pt != '<' )
2752                         ;
2753                 *pt = 0;
2754
2755
2756                 while( fgets( al, N-100, fp ) )
2757                         if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
2758
2759                 strcpy( tal, al+24 );
2760                 pt = tal;
2761                 while( *++pt != '<' )
2762                         ;
2763                 *pt = 0;
2764
2765
2766 //              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 );
2767
2768                 nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
2769
2770                 while( fgets( b, B-1, fp ) )
2771                         if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
2772
2773
2774                 fgets( b, B-1, fp );
2775
2776
2777                 if( !strncmp( "          </Hit_hsps>", b, 21 ) )
2778                 {
2779                         dis[junban[count++]] = sumscore;
2780                         sumscore = 0.0;
2781                         fgets( b, B-1, fp );
2782                         fgets( b, B-1, fp );
2783                         if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
2784                 }
2785         }
2786     return count;
2787 }
2788
2789 int ReadFasta34noalign( FILE *fp, double *dis, int qmem, char name[M][B], LocalHom *localhomlist )
2790 {
2791     int count=0;
2792     char b[B];
2793         char *pt;
2794     static int junban[M];
2795         int opt;
2796         double z, bits;
2797
2798
2799     count = 0;
2800 #if 0
2801     for( i=0; i<10000000 && count<nseq; i++ )
2802 #else
2803     while( !feof( fp ) )
2804 #endif
2805     {
2806         fgets( b, B-1, fp );
2807         if( !strncmp( "+==========+", b, 12 ) )
2808         {
2809             junban[count] = atoi( b+12 );
2810
2811                         pt = strchr( b, ')' ) + 1;
2812                         sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
2813             dis[junban[count]] = (double)opt;
2814             count++;
2815
2816         }
2817     }
2818
2819     return count;
2820 }
2821 int ReadFasta34m10_nuc( FILE *fp, double *dis, int qmem, char name[M][B], LocalHom *localhomlist )
2822 {
2823     int count=0;
2824     char b[B];
2825         char *pt;
2826     static int junban[M];
2827         int overlapaa;
2828         int opt, qstart, qend, tstart, tend;
2829         double z, bits;
2830         int qal_display_start, tal_display_start;
2831         static char qal[N], tal[N];
2832         char *qal2, *tal2;
2833         int c;
2834
2835
2836     count = 0;
2837 #if 0
2838     for( i=0; i<10000000 && count<nseq; i++ )
2839 #else
2840     while( !feof( fp ) )
2841 #endif
2842     {
2843         fgets( b, B-1, fp );
2844         if( !strncmp( "+==========+", b, 12 ) )
2845         {
2846             junban[count] = atoi( b+12 );
2847
2848                         if( strchr( b, 'r' ) ) continue;
2849
2850                         pt = strchr( b, ']' ) + 1;
2851                         sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
2852             dis[junban[count]] = (double)opt;
2853             count++;
2854
2855         }
2856                 else if( 0 == strncmp( ">>+==========+", b, 14 ) )
2857                 {
2858                         break;
2859                 }
2860
2861     }
2862         if( !count ) return -1;
2863
2864         count = 0;
2865     while( 1 )
2866         {
2867                 if( strncmp( ">>+==========+", b, 14 ) )
2868                 {
2869                         fgets( b, B-1, fp );
2870                         if( feof( fp ) ) break;
2871                         continue;
2872                 }
2873                 junban[count++] = atoi( b+14 );
2874 //              fprintf( stderr, "t = %d\n", atoi( b+14 ) );
2875                 while( fgets( b, B-1, fp ) )
2876                         if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
2877                 pt = strstr( b, ":" ) +1;
2878                 opt = atoi( pt );
2879
2880
2881                 while( fgets( b, B-1, fp ) )
2882                         if( !strncmp( "_overlap:", b+4, 9 ) ) break;
2883                 pt = strstr( b, ":" ) +1;
2884                 overlapaa = atoi( pt );
2885
2886                 while( fgets( b, B-1, fp ) )
2887                         if( !strncmp( "_start:", b+4, 7 ) ) break;
2888                 pt = strstr( b, ":" ) +1;
2889                 qstart = atoi( pt ) - 1;
2890
2891                 while( fgets( b, B-1, fp ) )
2892                         if( !strncmp( "_stop:", b+4, 6 ) ) break;
2893                 pt = strstr( b, ":" ) +1;
2894                 qend = atoi( pt ) - 1;
2895
2896                 while( fgets( b, B-1, fp ) )
2897                         if( !strncmp( "_display_start:", b+4, 15 ) ) break;
2898                 pt = strstr( b, ":" ) +1;
2899                 qal_display_start = atoi( pt ) - 1;
2900
2901                 pt = qal;
2902                 while( (c = fgetc( fp )) )
2903                 {
2904                         if( c == '>' )
2905                         {
2906                                 ungetc( c, fp );
2907                                 break;
2908                         }
2909                         if( isalpha( c ) || c == '-' )
2910                         *pt++ = c;
2911                 }
2912                 *pt = 0;
2913
2914                 while( fgets( b, B-1, fp ) )
2915                         if( !strncmp( "_start:", b+4, 7 ) ) break;
2916                 pt = strstr( b, ":" ) + 1;
2917                 tstart = atoi( pt ) - 1;
2918
2919                 while( fgets( b, B-1, fp ) )
2920                         if( !strncmp( "_stop:", b+4, 6 ) ) break;
2921                 pt = strstr( b, ":" ) + 1;
2922                 tend = atoi( pt ) - 1;
2923
2924                 while( fgets( b, B-1, fp ) )
2925                         if( !strncmp( "_display_start:", b+4, 15 ) ) break;
2926                 pt = strstr( b, ":" ) + 1;
2927                 tal_display_start = atoi( pt ) - 1;
2928
2929                 pt = tal;
2930                 while( ( c = fgetc( fp ) ) )
2931                 {
2932                         if( c == '>' )
2933                         {
2934                                 ungetc( c, fp );
2935                                 break;
2936                         }
2937                         if( isalpha( c ) || c == '-' )
2938                         *pt++ = c;
2939                 }
2940                 *pt = 0;
2941
2942 //              fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
2943 //              fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
2944
2945 //              fprintf( stderr, "qal = %s\n", qal );
2946 //              fprintf( stderr, "tal = %s\n", tal );
2947
2948                 qal2 = cutal( qal, qal_display_start, qstart, qend );
2949                 tal2 = cutal( tal, tal_display_start, tstart, tend );
2950
2951 //              fprintf( stderr, "qal2 = %s\n", qal2 );
2952 //              fprintf( stderr, "tal2 = %s\n", tal2 );
2953
2954 //              fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
2955                 putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
2956         }
2957 //      fprintf( stderr, "count = %d\n", count );
2958     return count;
2959 }
2960 int ReadFasta34m10( 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         int opt, qstart, qend, tstart, tend;
2968         double z, bits;
2969         int qal_display_start, tal_display_start;
2970         static char qal[N], tal[N];
2971         char *qal2, *tal2;
2972         int c;
2973
2974
2975     count = 0;
2976 #if 0
2977     for( i=0; i<10000000 && count<nseq; i++ )
2978 #else
2979     while( !feof( fp ) )
2980 #endif
2981     {
2982         fgets( b, B-1, fp );
2983         if( !strncmp( "+==========+", b, 12 ) )
2984         {
2985             junban[count] = atoi( b+12 );
2986
2987                         pt = strchr( b, ')' ) + 1;
2988                         sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
2989             dis[junban[count]] = (double)opt;
2990             count++;
2991
2992         }
2993                 else if( 0 == strncmp( ">>+==========+", b, 14 ) )
2994                 {
2995                         break;
2996                 }
2997
2998     }
2999         if( !count ) return -1;
3000
3001         count = 0;
3002     while( 1 )
3003         {
3004                 if( strncmp( ">>+==========+", b, 14 ) )
3005                 {
3006                         fgets( b, B-1, fp );
3007                         if( feof( fp ) ) break;
3008                         continue;
3009                 }
3010                 junban[count++] = atoi( b+14 );
3011 //              fprintf( stderr, "t = %d\n", atoi( b+14 ) );
3012                 while( fgets( b, B-1, fp ) )
3013                         if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
3014                 pt = strstr( b, ":" ) +1;
3015                 opt = atoi( pt );
3016
3017
3018                 while( fgets( b, B-1, fp ) )
3019                         if( !strncmp( "_overlap:", b+4, 9 ) ) break;
3020                 pt = strstr( b, ":" ) +1;
3021                 overlapaa = atoi( pt );
3022
3023                 while( fgets( b, B-1, fp ) )
3024                         if( !strncmp( "_start:", b+4, 7 ) ) break;
3025                 pt = strstr( b, ":" ) +1;
3026                 qstart = atoi( pt ) - 1;
3027
3028                 while( fgets( b, B-1, fp ) )
3029                         if( !strncmp( "_stop:", b+4, 6 ) ) break;
3030                 pt = strstr( b, ":" ) +1;
3031                 qend = atoi( pt ) - 1;
3032
3033                 while( fgets( b, B-1, fp ) )
3034                         if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3035                 pt = strstr( b, ":" ) +1;
3036                 qal_display_start = atoi( pt ) - 1;
3037
3038                 pt = qal;
3039                 while( (c = fgetc( fp )) )
3040                 {
3041                         if( c == '>' )
3042                         {
3043                                 ungetc( c, fp );
3044                                 break;
3045                         }
3046                         if( isalpha( c ) || c == '-' )
3047                         *pt++ = c;
3048                 }
3049                 *pt = 0;
3050
3051                 while( fgets( b, B-1, fp ) )
3052                         if( !strncmp( "_start:", b+4, 7 ) ) break;
3053                 pt = strstr( b, ":" ) + 1;
3054                 tstart = atoi( pt ) - 1;
3055
3056                 while( fgets( b, B-1, fp ) )
3057                         if( !strncmp( "_stop:", b+4, 6 ) ) break;
3058                 pt = strstr( b, ":" ) + 1;
3059                 tend = atoi( pt ) - 1;
3060
3061                 while( fgets( b, B-1, fp ) )
3062                         if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3063                 pt = strstr( b, ":" ) + 1;
3064                 tal_display_start = atoi( pt ) - 1;
3065
3066                 pt = tal;
3067                 while( ( c = fgetc( fp ) ) )
3068                 {
3069                         if( c == '>' )
3070                         {
3071                                 ungetc( c, fp );
3072                                 break;
3073                         }
3074                         if( isalpha( c ) || c == '-' )
3075                         *pt++ = c;
3076                 }
3077                 *pt = 0;
3078
3079 //              fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
3080 //              fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
3081
3082 //              fprintf( stderr, "qal = %s\n", qal );
3083 //              fprintf( stderr, "tal = %s\n", tal );
3084
3085                 qal2 = cutal( qal, qal_display_start, qstart, qend );
3086                 tal2 = cutal( tal, tal_display_start, tstart, tend );
3087
3088 //              fprintf( stderr, "qal2 = %s\n", qal2 );
3089 //              fprintf( stderr, "tal2 = %s\n", tal2 );
3090
3091 //              fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
3092                 putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
3093         }
3094 //      fprintf( stderr, "count = %d\n", count );
3095     return count;
3096 }
3097 int ReadFasta34m10_scoreonly_nucbk( FILE *fp, double *dis, int nin )
3098 {
3099     int count=0;
3100     char b[B];
3101         char *pt;
3102     int pos;
3103         int opt;
3104         double z, bits;
3105
3106     count = 0;
3107     while( !feof( fp ) )
3108     {
3109         fgets( b, B-1, fp );
3110         if( !strncmp( "+===========+", b, 13 ) )
3111         {
3112             pos = atoi( b+13 );
3113
3114                         if( strchr( b, 'r' ) ) continue;
3115
3116 //                      pt = strchr( b, ')' ) + 1;
3117                         pt = strchr( b, ']' ) + 1;
3118                         sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3119             dis[pos] += (double)opt;
3120             count++;
3121 #if 0
3122                         fprintf( stderr, "b=%s\n", b );
3123                         fprintf( stderr, "opt=%d\n", opt );
3124                         fprintf( stderr, "pos=%d\n", pos );
3125                         fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3126 #endif
3127
3128         }
3129                 else if( 0 == strncmp( ">>><<<", b, 6 ) )
3130                 {
3131                         break;
3132                 }
3133
3134     }
3135         if( !count ) return -1;
3136
3137     return count;
3138 }
3139
3140 int ReadFasta34m10_scoreonly_nuc( FILE *fp, double *dis, int nin )
3141 {
3142     int count=0;
3143     char b[B];
3144         char *pt;
3145     int pos;
3146         int opt;
3147         double z, bits;
3148         int c;
3149         int *yonda;
3150
3151
3152         yonda = AllocateIntVec( nin );
3153         for( c=0; c<nin; c++ ) yonda[c] = 0;
3154         for( c=0; c<nin; c++ ) dis[c] = 0.0;
3155
3156     count = 0;
3157     while( !feof( fp ) )
3158     {
3159         fgets( b, B-1, fp );
3160         if( !strncmp( "+===========+", b, 13 ) )
3161         {
3162             pos = atoi( b+13 );
3163
3164                         if( strchr( b, 'r' ) ) continue;
3165
3166 //                      pt = strchr( b, ')' ) + 1;
3167                         pt = strchr( b, ']' ) + 1;
3168                         sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3169                         if( yonda[pos] == 0 )
3170                         {
3171                     dis[pos] += (double)opt;
3172                                 yonda[pos] = 1;
3173                         }
3174             count++;
3175 #if 0
3176                         fprintf( stderr, "b=%s\n", b );
3177                         fprintf( stderr, "opt=%d\n", opt );
3178                         fprintf( stderr, "pos=%d\n", pos );
3179                         fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3180 #endif
3181
3182         }
3183         else if( !strncmp( ">>>", b, 3 ) )
3184                 {
3185                         for( c=0; c<nin; c++ ) yonda[c] = 0;
3186                 }
3187                 else if( 0 == strncmp( ">>><<<", b, 6 ) )
3188                 {
3189                         break;
3190                 }
3191
3192     }
3193
3194         free( yonda );
3195
3196         if( !count ) return -1;
3197
3198     return count;
3199 }
3200
3201 int ReadFasta34m10_scoreonly( FILE *fp, double *dis, int nin )
3202 {
3203     int count=0;
3204     char b[B];
3205         char *pt;
3206     int pos;
3207         int opt;
3208         double z, bits;
3209         int c;
3210         int *yonda;
3211
3212
3213         yonda = AllocateIntVec( nin );
3214         for( c=0; c<nin; c++ ) yonda[c] = 0;
3215         for( c=0; c<nin; c++ ) dis[c] = 0.0;
3216
3217     count = 0;
3218     while( !feof( fp ) )
3219     {
3220         fgets( b, B-1, fp );
3221         if( !strncmp( "+===========+", b, 13 ) )
3222         {
3223             pos = atoi( b+13 );
3224
3225                         pt = strchr( b, ')' ) + 1;
3226                         sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3227                         if( yonda[pos] == 0 )
3228                         {
3229                     dis[pos] += (double)opt;
3230                                 yonda[pos] = 1;
3231                         }
3232             count++;
3233 #if 0
3234                         fprintf( stderr, "b=%s\n", b );
3235                         fprintf( stderr, "opt=%d\n", opt );
3236                         fprintf( stderr, "pos=%d\n", pos );
3237                         fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3238 #endif
3239
3240         }
3241         else if( !strncmp( ">>>", b, 3 ) )
3242                 {
3243                         for( c=0; c<nin; c++ ) yonda[c] = 0;
3244                 }
3245                 else if( 0 == strncmp( ">>><<<", b, 6 ) )
3246                 {
3247                         break;
3248                 }
3249
3250     }
3251
3252         free( yonda );
3253
3254         if( !count ) return -1;
3255
3256     return count;
3257 }
3258 int ReadFasta34( FILE *fp, double *dis, int nseq, char name[M][B], LocalHom *localhomlist )
3259 {
3260     int count=0;
3261     char b[B];
3262         char *pt;
3263     static int junban[M];
3264         int overlapaa;
3265         int opt, qstart, qend, tstart, tend;
3266         double z, bits;
3267
3268
3269     count = 0;
3270 #if 0
3271     for( i=0; i<10000000 && count<nseq; i++ )
3272 #else
3273     while( !feof( fp ) )
3274 #endif
3275     {
3276         fgets( b, B-1, fp );
3277         if( !strncmp( "+==========+", b, 12 ) )
3278         {
3279             junban[count] = atoi( b+12 );
3280
3281                         pt = strchr( b, ')' ) + 1;
3282                         sscanf( pt, "%d %lf %lf",  &opt, &bits, &z );
3283             dis[junban[count]] = (double)opt;
3284             count++;
3285
3286         }
3287                 else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3288                 {
3289                         break;
3290                 }
3291
3292     }
3293         if( !count ) return -1;
3294
3295         count = 0;
3296     while( !feof( fp ) )
3297         {
3298                 if( !strncmp(">>+==========+", b, 14 ) )
3299                 {
3300             junban[count] = atoi( b+14 );
3301             count++;
3302                 fgets( b, B-1, fp ); // initn:
3303                         pt = strstr( b, "opt: " ) + 5;
3304                         localhomlist[junban[count-1]].opt = atof( pt );
3305                 fgets( b, B-1, fp ); // Smith-Waterman score
3306                         pt = strstr( b, "ungapped) in " ) + 13;
3307                         sscanf( pt, "%d", &overlapaa );
3308                         fprintf( stderr, "pt = %s, overlapaa = %d\n", pt, overlapaa );
3309                         pt = strstr( b, "overlap (" ) + 8;
3310                         sscanf( pt, "(%d-%d:%d-%d)", &qstart, &qend, &tstart, &tend );
3311                         localhomlist[junban[count-1]].overlapaa = overlapaa;
3312                         localhomlist[junban[count-1]].start1 = qstart-1;
3313                         localhomlist[junban[count-1]].end1   = qend-1;
3314                         localhomlist[junban[count-1]].start2 = tstart-1;
3315                         localhomlist[junban[count-1]].end2   = tend-1;
3316                 }
3317         fgets( b, B-1, fp );
3318         }
3319         fprintf( stderr, "count = %d\n", count );
3320     return count;
3321 }
3322
3323 int ReadFasta3( FILE *fp, double *dis, int nseq, char name[M][B] )
3324 {
3325     int count=0;
3326     char b[B];
3327         char *pt;
3328     int junban[M];
3329         int initn, init1, opt;
3330         double z;
3331
3332     count = 0;
3333 #if 0
3334     for( i=0; i<10000000 && count<nseq; i++ )
3335 #else
3336     while( !feof( fp ) )
3337 #endif
3338     {
3339         fgets( b, B-1, fp );
3340         if( !strncmp( "+==========+", b, 12 ) )
3341         {
3342             junban[count] = atoi( b+12 );
3343
3344                         pt = strchr( b, ')' ) + 1;
3345                         sscanf( pt, "%d %d %d %lf", &initn, &init1, &opt, &z );
3346             dis[junban[count]] = (double)opt;
3347             count++;
3348         }
3349     }
3350     return 0;
3351 }
3352
3353 int ReadFasta( FILE *fp, double *dis, int nseq, char name[M][B] )
3354 {
3355     int i, count=0;
3356     char b[B];
3357     int junban[M];
3358         int initn, init1, opt;
3359
3360     count = 0;
3361         for( i=0; i<nseq; i++ ) dis[i] = 0.0;
3362     for( i=0; !feof( fp ) && count<nseq; i++ )
3363     {
3364         fgets( b, B-1, fp );
3365         if( !strncmp( "+==========+", b, 12 ) )
3366         {
3367             junban[count] = atoi( b+12 );
3368
3369                         sscanf( b+50, "%d %d %d", &initn, &init1, &opt );
3370             dis[junban[count]] = (double)opt;
3371             count++;
3372         }
3373     }
3374
3375 /*
3376     count = 0;
3377     for( i=0; i<100000 && count<nseq; i++ )
3378     {
3379         fgets( b, B-1, fp );
3380         if( !strncmp( name[junban[count]], b, 20  ) )
3381         {
3382             dis[junban[count]] = atof( b+65 );
3383             count++;
3384         }
3385     }
3386 */
3387     return 0;
3388 }
3389
3390
3391 int ReadOpt( FILE *fp, int opt[M], int nseq, char name[M][B] )
3392 {
3393     int i, count=0;
3394     char b[B];
3395     int junban[M];
3396         int optt, initn, init1;
3397
3398     count = 0;
3399     for( i=0; i<10000000 && count<nseq; i++ )
3400     {
3401         fgets( b, B-1, fp );
3402         if( !strncmp( "+==========+", b, 12 ) )
3403         {
3404             junban[count] = atoi( b+12 );
3405                         sscanf( b+50, "%d %d %d", &initn, &init1, &optt );
3406             opt[junban[count]] = (double)optt;
3407             count++;
3408         }
3409     }
3410     return 0;
3411 }
3412
3413 int ReadOpt2( FILE *fp, int opt[M], int nseq, char name[M][B] )
3414 {
3415     int i, count=0;
3416     char b[B];
3417     int junban[M];
3418
3419     count = 0;
3420     for( i=0; i<10000000 && count<nseq; i++ )
3421     {
3422         fgets( b, B-1, fp );
3423         if( !strncmp( "+==========+", b, 12 ) )
3424         {
3425             junban[count] = atoi( b+12 );
3426             opt[junban[count]] = atoi( b+65 );
3427             count++;
3428         }
3429     }
3430     return 0;
3431 }
3432
3433
3434
3435 int writePre( int nseq, char name[][B], int nlen[M], char **aseq, int force )
3436 {
3437 #if USE_XCED
3438         int i, value;
3439         if( !signalSM )
3440         {
3441                 if( force )
3442                 {
3443                         rewind( prep_g );
3444                         if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3445 #if 0
3446                         else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3447 #else
3448                         else    writeData( prep_g, nseq, name, nlen, aseq );
3449 #endif
3450                 }
3451                 return( 0 );
3452         }
3453         for( i=0; i<10; i++ )
3454         {
3455 #if IODEBUG
3456                 fprintf( stderr, "SEMAPHORE = %d\n", signalSM[SEMAPHORE] );
3457 #endif
3458                 if( signalSM[SEMAPHORE]-- > 0 )
3459                 {
3460 #if 0 /* /tmp/pre ¤Î´Ø·¸¤Ç¤Ï¤º¤·¤¿ */
3461                         if( ferror( prep_g ) ) prep_g = fopen( "pre", "w" );
3462                         if( !prep_g ) ErrorExit( "Cannot re-open pre." );
3463 #endif
3464                         rewind( prep_g );
3465                         signalSM[STATUS] = IMA_KAITERU;
3466 #if IODEBUG
3467                         if( force ) fprintf( stderr, "FINAL " );
3468 #endif
3469                         if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3470                         else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3471                         /*
3472                         fprintf( prep_g, '\EOF' );
3473                         */
3474                         fflush( prep_g );
3475                         if( force ) signalSM[STATUS] = OSHIMAI;
3476                         else        signalSM[STATUS] = KAKIOWATTA;
3477                         value = 1;
3478                         signalSM[SEMAPHORE]++;
3479 #if IODEBUG
3480                         fprintf( stderr, "signalSM[STATUS] = %c\n", signalSM[STATUS] );
3481 #endif
3482                         break;
3483                 }
3484                 else
3485                 {
3486 #if IODEBUG
3487                         fprintf( stderr, "YONDERUKARA_AKIRAMERU\n" );
3488 #endif
3489                         value = 0;
3490                         signalSM[SEMAPHORE]++;
3491                         if( !force ) break;
3492 #if IODEBUG
3493                         fprintf( stderr, "MATSU\n" );
3494 #endif
3495                         sleep( 1 );
3496                 }
3497         }
3498         if( force && !value ) ErrorExit( "xced ga pre wo hanasanai \n" );
3499         return( value );
3500 #else
3501         if( force )
3502         {
3503                 rewind( prep_g );
3504                         writeData( prep_g, nseq, name, nlen, aseq );
3505         }
3506 #endif
3507         return( 0 );
3508 }
3509
3510
3511 void readOtherOptions( int *ppidptr, int *fftThresholdptr, int *fftWinSizeptr )
3512 {
3513         if( calledByXced )
3514         {
3515                 FILE *fp = fopen( "pre", "r" );
3516                 char b[B];
3517                 if( !fp ) ErrorExit( "Cannot open pre.\n" );
3518                 fgets( b, B-1, fp );
3519                 sscanf( b, "%d %d %d", ppidptr, fftThresholdptr, fftWinSizeptr );
3520                 fclose( fp );
3521 #if IODEBUG
3522         fprintf( stderr, "b = %s\n", b );
3523         fprintf( stderr, "ppid = %d\n", ppid );
3524         fprintf( stderr, "fftThreshold = %d\n", fftThreshold );
3525         fprintf( stderr, "fftWinSize = %d\n", fftWinSize );
3526 #endif
3527         }
3528         else
3529         {
3530                 *ppidptr = 0;
3531                 *fftThresholdptr = FFT_THRESHOLD;
3532                 if( dorp == 'd' )
3533                         *fftWinSizeptr = FFT_WINSIZE_D;
3534                 else
3535                         *fftWinSizeptr = FFT_WINSIZE_P;
3536         }
3537 #if 0
3538         fprintf( stderr, "fftThresholdptr=%d\n", *fftThresholdptr );
3539         fprintf( stderr, "fftWinSizeptr=%d\n", *fftWinSizeptr );
3540 #endif
3541 }
3542
3543 void initSignalSM( void )
3544 {
3545 //      int signalsmid;
3546
3547 #if IODEBUG
3548         if( ppid ) fprintf( stderr, "PID of xced = %d\n", ppid );
3549 #endif
3550         if( !ppid )
3551         {
3552                 signalSM = NULL;
3553                 return;
3554         }
3555
3556 #if 0
3557         signalsmid = shmget( (key_t)ppid, 3, IPC_ALLOC | 0666 );
3558         if( signalsmid == -1 ) ErrorExit( "Cannot get Shared memory for signal.\n" );
3559         signalSM = shmat( signalsmid, 0, 0 );
3560         if( (int)signalSM == -1 ) ErrorExit( "Cannot attatch Shared Memory for signal!\n" );
3561         signalSM[STATUS] = IMA_KAITERU;
3562         signalSM[SEMAPHORE] = 1;
3563 #endif
3564 }
3565
3566 void initFiles( void )
3567 {
3568         char pname[100];
3569         if( ppid )
3570                 sprintf( pname, "/tmp/pre.%d", ppid );
3571         else
3572                 sprintf( pname, "pre" );
3573         prep_g = fopen( pname, "w" );
3574         if( !prep_g ) ErrorExit( "Cannot open pre" );
3575
3576         trap_g = fopen( "trace", "w" );
3577         if( !trap_g ) ErrorExit( "cannot open trace" );
3578         fprintf( trap_g, "PID = %d\n", getpid() );
3579         fflush( trap_g );
3580 }
3581
3582
3583 void WriteForFasta( FILE *fp, int locnjob, char name[][B], int nlen[M], char **aseq )
3584 {
3585     static char b[N];
3586     int i, j;
3587     int nalen[M];
3588
3589     for( i=0; i<locnjob; i++ )
3590     {
3591         nalen[i] = strlen( aseq[i] );
3592         fprintf( fp, ">%s\n", name[i] );
3593         for( j=0; j<nalen[i]; j=j+C )
3594         {
3595             strncpy( b, aseq[i]+j, C ); b[C] = 0;
3596             fprintf( fp, "%s\n",b );
3597         }
3598     }
3599 }
3600
3601 void readlocalhomtable2( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
3602 {
3603         double opt;
3604         static char buff[B];
3605         char infor[100];
3606         int i, j, overlapaa, start1, end1, start2, end2;
3607         LocalHom *tmpptr1, *tmpptr2;
3608
3609 //      for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
3610
3611         while ( NULL != fgets( buff, B-1, fp ) )
3612         {
3613 //              fprintf( stderr, "\n" );
3614                 sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
3615                 if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
3616
3617 #if 0
3618                 if( start1 == end1 || start2 == end2 ) continue; //mondai ari
3619 #endif
3620
3621 //              if( i < j )
3622                 {
3623                         if( localhomtable[i][j].nokori++ > 0 )
3624                         {
3625                                 tmpptr1 = localhomtable[i][j].last;
3626 //                              fprintf( stderr, "reallocating, localhomtable[%d][%d].nokori = %d\n", i, j, localhomtable[i][j].nokori );
3627                                 tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
3628                                 tmpptr1 = tmpptr1->next;
3629                                 tmpptr1->extended = -1;
3630                                 tmpptr1->next = NULL;
3631                                 localhomtable[i][j].last = tmpptr1;
3632 //                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
3633                         }
3634                         else
3635                         {
3636                                 tmpptr1 = localhomtable[i]+j;
3637 //                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
3638                         }
3639
3640                         tmpptr1->start1 = start1;
3641                         tmpptr1->start2 = start2;
3642                         tmpptr1->end1 = end1;
3643                         tmpptr1->end2 = end2;
3644 //                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
3645 //                      tmpptr1->opt = opt;
3646                         tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
3647                         tmpptr1->overlapaa = overlapaa;
3648                         tmpptr1->korh = *infor;
3649                 }
3650 //              else
3651                 {
3652                         if( localhomtable[j][i].nokori++ > 0 )
3653                         {
3654                                 tmpptr2 = localhomtable[j][i].last;
3655                                 tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
3656                                 tmpptr2 = tmpptr2->next;
3657                                 tmpptr2->extended = -1;
3658                                 tmpptr2->next = NULL;
3659                                 localhomtable[j][i].last = tmpptr2;
3660 //                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
3661                         }
3662                         else
3663                         {
3664                                 tmpptr2 = localhomtable[j]+i;
3665 //                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
3666                         }
3667
3668                         tmpptr2->start2 = start1;
3669                         tmpptr2->start1 = start2;
3670                         tmpptr2->end2 = end1;
3671                         tmpptr2->end1 = end2;
3672 //                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
3673 //                      tmpptr2->opt = opt;
3674                         tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
3675                         tmpptr2->overlapaa = overlapaa;
3676                         tmpptr2->korh = *infor;
3677
3678 //                      fprintf( stderr, "i=%d, j=%d, st1=%d, en1=%d, opt = %f\n", i, j, tmpptr1->start1, tmpptr1->end1, opt );
3679                 }
3680
3681         }
3682 }
3683 void readlocalhomtable( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
3684 {
3685         double opt;
3686         static char buff[B];
3687         char infor[100];
3688         int i, j, overlapaa, start1, end1, start2, end2;
3689         int **nlocalhom = NULL;
3690         LocalHom *tmpptr1=NULL, *tmpptr2=NULL; // by D.Mathog, a guess
3691
3692         nlocalhom = AllocateIntMtx( njob, njob );
3693         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
3694
3695         while ( NULL != fgets( buff, B-1, fp ) )
3696         {
3697 //              fprintf( stderr, "\n" );
3698                 sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
3699                 if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
3700
3701 #if 0
3702                 if( start1 == end1 || start2 == end2 ) continue; //mondai ari
3703 #endif
3704
3705
3706 //              if( i < j )
3707                 {
3708                         if( nlocalhom[i][j]++ > 0 )
3709                         {
3710 //                              fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
3711                                 tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
3712                                 tmpptr1 = tmpptr1->next;
3713                                 tmpptr1->next = NULL;
3714                         }
3715                         else
3716                         {
3717                                 tmpptr1 = localhomtable[i]+j;
3718 //                              fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
3719                         }
3720
3721                         tmpptr1->start1 = start1;
3722                         tmpptr1->start2 = start2;
3723                         tmpptr1->end1 = end1;
3724                         tmpptr1->end2 = end2;
3725 //                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
3726 //                      tmpptr1->opt = opt;
3727                         tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
3728                         tmpptr1->overlapaa = overlapaa;
3729                         tmpptr1->korh = *infor;
3730
3731 //                      fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
3732                 }
3733 //              else
3734                 {
3735                         if( nlocalhom[j][i]++ > 0 )
3736                         {
3737                                 tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
3738                                 tmpptr2 = tmpptr2->next;
3739                                 tmpptr2->next = NULL;
3740                         }
3741                         else
3742                                 tmpptr2 = localhomtable[j]+i;
3743
3744                         tmpptr2->start2 = start1;
3745                         tmpptr2->start1 = start2;
3746                         tmpptr2->end2 = end1;
3747                         tmpptr2->end1 = end2;
3748 //                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
3749 //                      tmpptr2->opt = opt;
3750                         tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
3751                         tmpptr2->overlapaa = overlapaa;
3752                         tmpptr2->korh = *infor;
3753                 }
3754
3755         }
3756         FreeIntMtx( nlocalhom );
3757 }
3758
3759 void outlocalhom( LocalHom **localhom, int nseq )
3760 {
3761         int i, j;
3762         LocalHom *tmpptr;
3763         for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
3764         {
3765                 tmpptr = localhom[i]+j;
3766                 fprintf( stderr, "%d-%d\n", i, j );
3767                 do
3768                 {
3769                         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 );
3770                 }
3771                 while( (tmpptr=tmpptr->next) );
3772         }
3773 }
3774
3775 void outlocalhompt( LocalHom ***localhom, int n1, int n2 )
3776 {
3777         int i, j;
3778         LocalHom *tmpptr;
3779         for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
3780         {
3781                 tmpptr = localhom[i][j];
3782                 fprintf( stderr, "%d-%d\n", i, j );
3783                 do
3784                 {
3785                         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 );
3786                 }
3787                 while( (tmpptr=tmpptr->next) );
3788         }
3789 }
3790
3791 void FreeLocalHomTable( LocalHom **localhomtable, int n )
3792 {
3793         int i, j;
3794         LocalHom *ppp, *tmpptr;
3795         fprintf( stderr, "freeing localhom\n" );
3796         for( i=0; i<n; i++ )
3797         {
3798                 for( j=0; j<n; j++ )
3799                 {
3800                         tmpptr=localhomtable[i]+j;
3801                         ppp = tmpptr->next;
3802                         for( ; tmpptr; tmpptr=ppp )
3803                         {
3804 #if DEBUG
3805                                 fprintf( stderr, "i=%d, j=%d\n", i, j );
3806 #endif
3807                                 ppp = tmpptr->next;
3808                                 if( tmpptr!=localhomtable[i]+j )
3809                                 {
3810 #if DEBUG
3811                                         fprintf( stderr, "freeing %p\n", tmpptr );
3812 #endif
3813                                         free( tmpptr );
3814                                 }
3815                         }
3816                 }
3817 #if DEBUG
3818                 fprintf( stderr, "freeing localhomtable[%d]\n", i );
3819 #endif
3820                 free( localhomtable[i] );
3821         }
3822 #if DEBUG
3823         fprintf( stderr, "freeing localhomtable\n" );
3824 #endif
3825         free( localhomtable );
3826 #if DEBUG
3827         fprintf( stderr, "freed\n" );
3828 #endif
3829 }
3830
3831 char *progName( char *str )
3832 {
3833     char *value;
3834     if( ( value = strrchr( str, '/' ) ) != NULL )
3835         return( value+1 );
3836     else
3837         return( str );
3838 }
3839
3840 static void tabtospace( char *str )
3841 {
3842         char *p;
3843 //      fprintf( stderr, "before = %s\n", str );
3844         while( NULL != ( p = strchr( str , '\t' ) ) )
3845         {
3846                 *p = ' ';
3847         }
3848 //      fprintf( stderr, "after = %s\n", str );
3849 }
3850
3851 static char *extractfirstword( char *str )
3852 {
3853         char *val = str;
3854
3855         tabtospace( str );
3856         while( *str )
3857         {
3858                 if( val == str && *str == ' ' )
3859                 {
3860                         val++; str++;
3861                 }
3862                 else if( *str != ' ' )
3863                 {
3864                         str++;
3865                 }
3866                 else if( *str == ' ' )
3867                 {
3868                         *str = 0;
3869                 }
3870         }
3871         return( val );
3872 }
3873
3874
3875 void clustalout_pointer( FILE *fp, int nseq, int maxlen, char **seq, char **name, char *mark, char *comment, int *order )
3876 {
3877         int pos, j;
3878         pos = 0;
3879         if( comment == NULL )
3880                 fprintf( fp, "CLUSTAL format alignment by MAFFT (v%s)\n\n", VERSION );
3881         else
3882                 fprintf( fp, "CLUSTAL format alignment by MAFFT %s (v%s)\n\n", comment, VERSION );
3883
3884         while( pos < maxlen )
3885         {
3886                 fprintf( fp, "\n" );
3887                 for( j=0; j<nseq; j++ )
3888                 {
3889                         fprintf( fp, "%-15.15s ", extractfirstword( name[order[j]]+1 ) );
3890                         fprintf( fp, "%.60s\n", seq[order[j]]+pos ); // Ä¹¤µ¤¬°ã¤¦¤È¤À¤á¡£
3891                 }
3892                 if( mark )
3893                 {
3894                         fprintf( fp, "%-15.15s ", "" );
3895                         fprintf( fp, "%.60s\n", mark + pos ); // Ä¹¤µ¤¬°ã¤¦¤È¤À¤á¡£
3896                 }
3897                 pos += 60;
3898         }
3899 }
3900
3901
3902 void writeData_reorder_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq, int *order )
3903 {
3904         int i, j, k;
3905         int nalen;
3906
3907         for( i=0; i<locnjob; i++ )
3908         {
3909                 k = order[i];
3910 #if DEBUG
3911                 fprintf( stderr, "i = %d in writeData\n", i );
3912 #endif
3913                 nalen = strlen( aseq[k] );
3914                 fprintf( fp, ">%s\n", name[k]+1 );
3915                 for( j=0; j<nalen; j=j+C )
3916                 {
3917 #if 0
3918                         strncpy( b, aseq[k]+j, C ); b[C] = 0;
3919                         fprintf( fp, "%s\n",b );
3920 #else
3921                         fprintf( fp, "%.*s\n", C, aseq[k]+j );
3922 #endif
3923                 }
3924         }
3925 }
3926 void writeData_reorder( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq, int *order )
3927 {
3928         int i, j, k;
3929         int nalen;
3930
3931         for( i=0; i<locnjob; i++ )
3932         {
3933                 k = order[i];
3934 #if DEBUG
3935                 fprintf( stderr, "i = %d in writeData\n", i );
3936 #endif
3937                 nalen = strlen( aseq[k] );
3938                 fprintf( fp, ">%s\n", name[k]+1 );
3939                 for( j=0; j<nalen; j=j+C )
3940                 {
3941 #if 0
3942                         strncpy( b, aseq[k]+j, C ); b[C] = 0;
3943                         fprintf( fp, "%s\n",b );
3944 #else
3945                         fprintf( fp, "%.*s\n", C, aseq[k]+j );
3946 #endif
3947                 }
3948         }
3949 }
3950 static void showaamtxexample()
3951 {
3952         fprintf( stderr, "Format error in aa matrix\n" );
3953         fprintf( stderr, "# Example:\n" );
3954         fprintf( stderr, "# comment\n" );
3955         fprintf( stderr, "   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V\n" );
3956         fprintf( stderr, "A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0\n" );
3957         fprintf( stderr, "R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3\n" );
3958         fprintf( stderr, "...\n" );
3959         fprintf( stderr, "V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4\n" );
3960         fprintf( stderr, "frequency 0.07 0.05 0.04 0.05 0.02 .. \n" );
3961         fprintf( stderr, "# Example end\n" );
3962         fprintf( stderr, "Only the lower half is loaded\n" );
3963         fprintf( stderr, "The last line (frequency) is optional.\n" );
3964         exit( 1 );
3965 }
3966
3967 double *loadaamtx( void )
3968 {
3969         int i, j, k, ii, jj;
3970         double *val;
3971         double **raw;
3972         int *map;
3973         char *aaorder = "ARNDCQEGHILKMFPSTWYV";
3974         char *inorder;
3975         char *line;
3976         char *ptr1;
3977         char *ptr2;
3978         char *mtxfname = "_aamtx";
3979         FILE *mf;
3980
3981         raw = AllocateDoubleMtx( 21, 20 );
3982         val = AllocateDoubleVec( 420 );
3983         map = AllocateIntVec( 20 );
3984
3985         if( dorp != 'p' )
3986         {
3987                 fprintf( stderr, "User-defined matrix is not supported for DNA\n" );
3988                 exit( 1 );
3989         }
3990
3991         mf = fopen( mtxfname, "r" );
3992         if( mf == NULL )
3993         {
3994                 fprintf( stderr, "Cannot open the _aamtx file\n" );
3995                 exit( 1 );
3996         }
3997
3998         inorder = calloc( 1000, sizeof( char ) );
3999         line = calloc( 1000, sizeof( char ) );
4000
4001
4002         while( !feof( mf ) )
4003         {
4004                 fgets( inorder, 999, mf );
4005                 if( inorder[0] != '#' ) break;
4006         }
4007         ptr1 = ptr2 = inorder;
4008         while( *ptr2 )
4009         {
4010                 if( isalpha( *ptr2 ) )
4011                 {
4012                         *ptr1 = toupper( *ptr2 );
4013                         ptr1++;
4014                 }
4015                 ptr2++;
4016         }
4017         inorder[20] = 0;
4018
4019         for( i=0; i<20; i++ )
4020         {
4021                 ptr2 = strchr( inorder, aaorder[i] );
4022                 if( ptr2 == NULL )
4023                 {
4024                         fprintf( stderr, "%c: not found in the first 20 letters.\n", aaorder[i] );
4025                         showaamtxexample();
4026                 }
4027                 else
4028                 {
4029                         map[i] = ptr2 - inorder;
4030                 }
4031         }
4032
4033         i = 0;
4034         while( !feof( mf ) )
4035         {
4036                 fgets( line, 999, mf );
4037 //              fprintf( stderr, "line = %s\n", line );
4038                 if( line[0] == '#' ) continue;
4039                 ptr1 = line;
4040 //              fprintf( stderr, "line = %s\n", line );
4041                 for( j=0; j<=i; j++ )
4042                 {
4043                         while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4044                                 ptr1++;
4045
4046                         raw[i][j] = atof( ptr1 );
4047 //                      fprintf( stderr, "raw[][]=%f, %c-%c %d-%d\n", raw[i][j], inorder[i], inorder[j], i, j );
4048                         ptr1 = strchr( ptr1, ' ' );
4049                         if( ptr1 == NULL && j<i) showaamtxexample();
4050                 }
4051                 i++;
4052                 if( i > 19 ) break;
4053         }
4054
4055         for( i=0; i<20; i++ ) raw[20][i] = -1.0;
4056         while( !feof( mf ) )
4057         {
4058                 fgets( line, 999, mf );
4059                 if( line[0] == 'f' )
4060                 {
4061 //                      fprintf( stderr, "line = %s\n", line );
4062                         ptr1 = line;
4063                         for( j=0; j<20; j++ )
4064                         {
4065                                 while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4066                                         ptr1++;
4067
4068                                 raw[20][j] = atof( ptr1 );
4069 //                              fprintf( stderr, "raw[20][]=%f, %c %d\n", raw[20][j], inorder[i], j );
4070                                 ptr1 = strchr( ptr1, ' ' );
4071                                 if( ptr1 == NULL && j<19) showaamtxexample();
4072                         }
4073                         break;
4074                 }
4075         }
4076
4077         k = 0;
4078         for( i=0; i<20; i++ )
4079         {
4080                 for( j=0; j<=i; j++ )
4081                 {
4082                         if( i != j )
4083                         {
4084                                 ii = MAX( map[i], map[j] );
4085                                 jj = MIN( map[i], map[j] );
4086                         }
4087                         else ii = jj = map[i];
4088                         val[k++] = raw[ii][jj];
4089 //                      fprintf( stderr, "%c-%c, %f\n", aaorder[i], aaorder[j], val[k-1] );
4090                 }
4091         }
4092         for( i=0; i<20; i++ ) val[400+i] = raw[20][map[i]];
4093
4094         fprintf( stderr, "inorder = %s\n", inorder );
4095         fclose( mf );
4096         free( inorder );
4097         free( line );
4098         FreeDoubleMtx( raw );
4099         free( map );
4100         return( val );
4101 }
4102
4103 void miyataout_reorder_pointer( FILE *fp, int locnjob, int nlenmax, char **name, int *nlen, char **aseq, int *order )
4104 {
4105         int i, j, k;
4106         int nalen;
4107
4108         fprintf( fp, "%5d\n", 1 );
4109         fprintf( fp, "%5d\n", 1 );
4110         fprintf( fp, "%5d%5d\n", 1, nlenmax );
4111         fprintf( fp, "%5d\n", 0 );
4112         fprintf( fp, "%5d\n", locnjob );
4113         for( i=0; i<locnjob; i++ )
4114         {
4115                 k = order[i];
4116                 nalen = strlen( aseq[k] );
4117                 fprintf( fp, "=%s\n%d\n", name[k]+1, nalen );
4118                 for( j=0; j<nalen; j=j+C )
4119                 {
4120                         fprintf( fp, "%.*s\n", C, aseq[k]+j );
4121                 }
4122         }
4123 }
4124
4125 void readmccaskill( FILE *fp, RNApair **pairprob, int length )
4126 {
4127         char gett[1000];
4128         int *pairnum;
4129         int i;
4130         int left, right;
4131         float prob;
4132         int c;
4133
4134         pairnum = (int *)calloc( length, sizeof( int ) );
4135         for( i=0; i<length; i++ ) pairnum[i] = 0;
4136
4137         c = getc( fp );
4138         {
4139                 if( c != '>' )
4140                 {
4141                         fprintf( stderr, "format error in hat4\n" );
4142                         exit( 1 );
4143                 }
4144         }
4145         fgets( gett, 999, fp );
4146         while( 1 )
4147         {
4148                 if( feof( fp ) ) break;
4149                 c = getc( fp );
4150                 ungetc( c, fp );
4151                 if( c == '>' )
4152                 {
4153                         break;
4154                 }
4155                 fgets( gett, 999, fp );
4156 //              fprintf( stderr, "gett = %s\n", gett );
4157                 sscanf( gett, "%d %d %f", &left, &right, &prob );
4158
4159                 if( left >= length || right >= length )
4160                 {
4161                         fprintf( stderr, "format error in hat4\n" );
4162                         exit( 1 );
4163                 }
4164
4165                 if( prob < 0.01 ) continue; // 080607, mafft ni dake eikyou
4166
4167                 if( left != right && prob > 0.0 )
4168                 {
4169                         pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
4170                         pairprob[left][pairnum[left]].bestscore = prob;
4171                         pairprob[left][pairnum[left]].bestpos = right;
4172                         pairnum[left]++;
4173                         pairprob[left][pairnum[left]].bestscore = -1.0;
4174                         pairprob[left][pairnum[left]].bestpos = -1;
4175 //                      fprintf( stderr, "%d-%d, %f\n", left, right, prob );
4176
4177                         pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
4178                         pairprob[right][pairnum[right]].bestscore = prob;
4179                         pairprob[right][pairnum[right]].bestpos = left;
4180                         pairnum[right]++;
4181                         pairprob[right][pairnum[right]].bestscore = -1.0;
4182                         pairprob[right][pairnum[right]].bestpos = -1;
4183 //                      fprintf( stderr, "%d-%d, %f\n", right, left, prob );
4184                 }
4185         }
4186         free( pairnum );
4187 }
4188
4189 void readpairfoldalign( FILE *fp, char *s1, char *s2, char *aln1, char *aln2, int q1, int q2, int *of1, int *of2, int sumlen )
4190 {
4191         char gett[1000];
4192         int *maptoseq1;
4193         int *maptoseq2;
4194         char dumc;
4195         int dumi;
4196         char sinseq[100], sinaln[100];
4197         int posinseq, posinaln;
4198         int alnlen;
4199         int i;
4200         int pos1, pos2;
4201         char *pa1, *pa2;
4202         char qstr[1000];
4203
4204         *of1 = -1;
4205         *of2 = -1;
4206
4207         maptoseq1 = AllocateIntVec( sumlen+1 );
4208         maptoseq2 = AllocateIntVec( sumlen+1 );
4209
4210         posinaln = 0; // foldalign ga alingment wo kaesanaitok no tame.
4211
4212         while( !feof( fp ) )
4213         {
4214                 fgets( gett, 999, fp );
4215                 if( !strncmp( gett, "; ALIGNING", 10 ) ) break;
4216         }
4217         sprintf( qstr, "; ALIGNING            %d against %d\n", q1+1, q2+1 );
4218         if( strcmp( gett, qstr ) )
4219         {
4220                 fprintf( stderr, "Error in FOLDALIGN\n" );
4221                 fprintf( stderr, "qstr = %s, but gett = %s\n", qstr, gett );
4222                 exit( 1 );
4223         }
4224
4225         while( !feof( fp ) )
4226         {
4227                 fgets( gett, 999, fp );
4228                 if( !strncmp( gett, "; --------", 10 ) ) break;
4229         }
4230
4231
4232         while( !feof( fp ) )
4233         {
4234                 fgets( gett, 999, fp );
4235                 if( !strncmp( gett, "; ********", 10 ) ) break;
4236 //              fprintf( stderr, "gett = %s\n", gett );
4237                 sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
4238                 posinaln = atoi( sinaln );
4239                 posinseq = atoi( sinseq );
4240 //              fprintf( stderr, "posinseq = %d\n", posinseq );
4241 //              fprintf( stderr, "posinaln = %d\n", posinaln );
4242                 maptoseq1[posinaln-1] = posinseq-1;
4243         }
4244         alnlen = posinaln;
4245
4246         while( !feof( fp ) )
4247         {
4248                 fgets( gett, 999, fp );
4249                 if( !strncmp( gett, "; --------", 10 ) ) break;
4250         }
4251
4252         while( !feof( fp ) )
4253         {
4254                 fgets( gett, 999, fp );
4255                 if( !strncmp( gett, "; ********", 10 ) ) break;
4256 //              fprintf( stderr, "gett = %s\n", gett );
4257                 sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
4258                 posinaln = atof( sinaln );
4259                 posinseq = atof( sinseq );
4260 //              fprintf( stderr, "posinseq = %d\n", posinseq );
4261 //              fprintf( stderr, "posinaln = %d\n", posinaln );
4262                 maptoseq2[posinaln-1] = posinseq-1;
4263         }
4264         if( alnlen != posinaln )
4265         {
4266                 fprintf( stderr, "Error in foldalign?\n" );
4267                 exit( 1 );
4268         }
4269
4270         pa1 = aln1;
4271         pa2 = aln2;
4272         for( i=0; i<alnlen; i++ )
4273         {
4274                 pos1 = maptoseq1[i];
4275                 pos2 = maptoseq2[i];
4276
4277                 if( pos1 > -1 )
4278                         *pa1++ = s1[pos1];
4279                 else
4280                         *pa1++ = '-';
4281
4282                 if( pos2 > -1 )
4283                         *pa2++ = s2[pos2];
4284                 else
4285                         *pa2++ = '-';
4286         }
4287         *pa1 = 0;
4288         *pa2 = 0;
4289
4290         *of1 = 0;
4291         for( i=0; i<alnlen; i++ )
4292         {
4293                 *of1 = maptoseq1[i];
4294                 if( *of1 > -1 ) break;
4295         }
4296         *of2 = 0;
4297         for( i=0; i<alnlen; i++ )
4298         {
4299                 *of2 = maptoseq2[i];
4300                 if( *of2 > -1 ) break;
4301         }
4302
4303 //      fprintf( stderr, "*of1=%d, aln1 = :%s:\n", *of1, aln1 );
4304 //      fprintf( stderr, "*of2=%d, aln2 = :%s:\n", *of2, aln2 );
4305
4306         free( maptoseq1 );
4307         free( maptoseq2 );
4308 }