binaries/src/mafft/core/tbfast3.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 0
   6
   7
   8 void arguments( int argc, char *argv[] )
   9 {
  10     int c;
  11
  12         thrinter = 1.0;
  13         inputfile = NULL;
  14         fftkeika = 0;
  15         constraint = 0;
  16         nblosum = 62;
  17         fmodel = 0;
  18         calledByXced = 0;
  19         devide = 0;
  20         use_fft = 0;
  21         fftscore = 1;
  22         fftRepeatStop = 0;
  23         fftNoAnchStop = 0;
  24     weight = 3;
  25     utree = 1;
  26         tbutree = 1;
  27     refine = 0;
  28     check = 1;
  29     cut = 0.0;
  30     disp = 0;
  31     outgap = 1;
  32     alg = 'A';
  33     mix = 0;
  34         tbitr = 0;
  35         scmtd = 5;
  36         tbweight = 0;
  37         tbrweight = 3;
  38         checkC = 0;
  39         treemethod = 'x';
  40         contin = 0;
  41         scoremtx = 1;
  42         kobetsubunkatsu = 0;
  43         dorp = NOTSPECIFIED;
  44         ppenalty = NOTSPECIFIED;
  45         ppenalty_ex = NOTSPECIFIED;
  46         poffset = NOTSPECIFIED;
  47         kimuraR = NOTSPECIFIED;
  48         pamN = NOTSPECIFIED;
  49         geta2 = GETA2;
  50         fftWinSize = NOTSPECIFIED;
  51         fftThreshold = NOTSPECIFIED;
  52         TMorJTT = JTT;
  53
  54     while( --argc > 0 && (*++argv)[0] == '-' )
  55         {
  56         while ( c = *++argv[0] )
  57                 {
  58             switch( c )
  59             {
  60                                 case 'i':
  61                                         inputfile = *++argv;
  62                                         fprintf( stderr, "inputfile = %s\n", inputfile );
  63                                         --argc;
  64                     goto nextoption;
  65                                 case 'f':
  66                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
  67                                         fprintf( stderr, "ppenalty = %d\n", ppenalty );
  68                                         --argc;
  69                                         goto nextoption;
  70                                 case 'g':
  71                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
  72                                         fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
  73                                         --argc;
  74                                         goto nextoption;
  75                                 case 'h':
  76                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
  77                                         fprintf( stderr, "poffset = %d\n", poffset );
  78                                         --argc;
  79                                         goto nextoption;
  80                                 case 'k':
  81                                         kimuraR = atoi( *++argv );
  82                                         fprintf( stderr, "kimuraR = %d\n", kimuraR );
  83                                         --argc;
  84                                         goto nextoption;
  85                                 case 'b':
  86                                         nblosum = atoi( *++argv );
  87                                         scoremtx = 1;
  88                                         fprintf( stderr, "blosum %d\n", nblosum );
  89                                         --argc;
  90                                         goto nextoption;
  91                                 case 'j':
  92                                         pamN = atoi( *++argv );
  93                                         scoremtx = 0;
  94                                         TMorJTT = JTT;
  95                                         fprintf( stderr, "jtt %d\n", pamN );
  96                                         --argc;
  97                                         goto nextoption;
  98                                 case 'm':
  99                                         pamN = atoi( *++argv );
 100                                         scoremtx = 0;
 101                                         TMorJTT = TM;
 102                                         fprintf( stderr, "tm %d\n", pamN );
 103                                         --argc;
 104                                         goto nextoption;
 105                                 case 'l':
 106                                         fastathreshold = atof( *++argv );
 107                                         constraint = 2;
 108                                         fprintf( stderr, "weighti = %f\n", fastathreshold );
 109                                         --argc;
 110                                         goto nextoption;
 111                                 case 'c':
 112                                         thrinter = atof( *++argv );
 113                                         fprintf( stderr, "thrinter %f\n", thrinter );
 114                                         --argc;
 115                                         goto nextoption;
 116 #if 0
 117                                 case 'm':
 118                                         fmodel = 1;
 119                                         break;
 120 #endif
 121                                 case 'r':
 122                                         fmodel = -1;
 123                                         break;
 124                                 case 'D':
 125                                         dorp = 'd';
 126                                         break;
 127                                 case 'P':
 128                                         dorp = 'p';
 129                                         break;
 130                                 case 'e':
 131                                         fftscore = 0;
 132                                         break;
 133                                 case 'O':
 134                                         fftNoAnchStop = 1;
 135                                         break;
 136                                 case 'R':
 137                                         fftRepeatStop = 1;
 138                                         break;
 139                                 case 'Q':
 140                                         calledByXced = 1;
 141                                         break;
 142                                 case 's':
 143                                         treemethod = 's';
 144                                         break;
 145                                 case 'x':
 146                                         treemethod = 'x';
 147                                         break;
 148                                 case 'p':
 149                                         treemethod = 'p';
 150                                         break;
 151                                 case 'a':
 152                                         alg = 'a';
 153                                         break;
 154                                 case 'A':
 155                                         alg = 'A';
 156                                         break;
 157                                 case 'S':
 158                                         alg = 'S';
 159                                         break;
 160                                 case 'C':
 161                                         alg = 'C';
 162                                         break;
 163                                 case 'F':
 164                                         use_fft = 1;
 165                                         break;
 166                                 case 'v':
 167                                         tbrweight = 3;
 168                                         break;
 169                                 case 'd':
 170                                         disp = 1;
 171                                         break;
 172                                 case 'o':
 173                                         outgap = 0;
 174                                         break;
 175 /* Modified 01/08/27, default: user tree */
 176                                 case 'J':
 177                                         tbutree = 0;
 178                                         break;
 179 /* modification end. */
 180                                 case 'z':
 181                                         fftThreshold = atoi( *++argv );
 182                                         --argc;
 183                                         goto nextoption;
 184                                 case 'w':
 185                                         fftWinSize = atoi( *++argv );
 186                                         --argc;
 187                                         goto nextoption;
 188                                 case 'Z':
 189                                         checkC = 1;
 190                                         break;
 191                 default:
 192                     fprintf( stderr, "illegal option %c\n", c );
 193                     argc = 0;
 194                     break;
 195             }
 196                 }
 197                 nextoption:
 198                         ;
 199         }
 200     if( argc == 1 )
 201     {
 202         cut = atof( (*argv) );
 203         argc--;
 204     }
 205     if( argc != 0 )
 206     {
 207         fprintf( stderr, "options: Check source file !\n" );
 208         exit( 1 );
 209     }
 210         if( tbitr == 1 && outgap == 0 )
 211         {
 212                 fprintf( stderr, "conflicting options : o, m or u\n" );
 213                 exit( 1 );
 214         }
 215         if( alg == 'C' && outgap == 0 )
 216         {
 217                 fprintf( stderr, "conflicting options : C, o\n" );
 218                 exit( 1 );
 219         }
 220 }
 221
 222
 223 void treebase( char name[M][B], int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, int ***topol, double *effarr, int alloclen, LocalHom **localhomtable )
 224 {
 225         int i, j, l;
 226         int clus1, clus2;
 227         int s1, s2, r1, r2;
 228         float pscore, tscore;
 229         static char *indication1, *indication2;
 230         FILE *trap;
 231         static char name1[M][B], name2[M][B];
 232         static double *effarr1 = NULL;
 233         static double *effarr2 = NULL;
 234         static LocalHom ***localhomshrink = NULL;
 235         float dumfl = 0.0;
 236         double total;
 237
 238
 239         if( effarr1 == NULL )
 240         {
 241                 effarr1 = AllocateDoubleVec( njob );
 242                 effarr2 = AllocateDoubleVec( njob );
 243                 indication1 = AllocateCharVec( 150 );
 244                 indication2 = AllocateCharVec( 150 );
 245 #if 0
 246 #else
 247                 if( constraint )
 248                 {
 249                         localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
 250                         for( i=0; i<njob; i++)
 251                                 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 252                 }
 253 #endif
 254         }
 255
 256 #if 0
 257         fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
 258 #endif
 259
 260 #if 0
 261         for( i=0; i<njob; i++ )
 262                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
 263 #endif
 264
 265         for( i=0; i<njob; i++ ) strcpy( aseq[i], seq[i] );
 266
 267         if( constraint )
 268                 calcimportance( njob, effarr, aseq, localhomtable );
 269
 270
 271         writePre( njob, name, nlen, aseq, 0 );
 272
 273         tscore = 0.0;
 274         for( l=0; l<njob-1; l++ )
 275         {
 276
 277                 clus1 = fastconjuction( topol[l][0], aseq, mseq1, effarr1, effarr, name, name1, indication1 );
 278                 clus2 = fastconjuction( topol[l][1], aseq, mseq2, effarr2, effarr, name, name2, indication2 );
 279
 280
 281                 fprintf( trap_g, "\nSTEP-%d\n", l );
 282                 fprintf( trap_g, "group1 = %s\n", indication1 );
 283                 fprintf( trap_g, "group2 = %s\n", indication2 );
 284
 285 #if 1
 286                 fprintf( stderr, "STEP %d /%d\r", l+1, njob-1 );
 287 #else
 288                 fprintf( stderr, "STEP %d /%d\n", l+1, njob-1 );
 289                 fprintf( stderr, "group1 = %.66s", indication1 );
 290                 if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
 291                 fprintf( stderr, "\n" );
 292                 fprintf( stderr, "group2 = %.66s", indication2 );
 293                 if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
 294                 fprintf( stderr, "\n" );
 295 #endif
 296 //              for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
 297
 298                 if( constraint )
 299                 {
 300                         fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
 301 //                      fprintf( stderr, "localhomshrink =\n" );
 302 //                      outlocalhompt( localhomshrink, clus1, clus2 );
 303 //                      weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
 304 //                      fprintf( stderr, "after weight =\n" );
 305 //                      outlocalhompt( localhomshrink, clus1, clus2 );
 306                 }
 307
 308 /*
 309                 fprintf( stderr, "before align all\n" );
 310                 display( aseq, njob );
 311                 fprintf( stderr, "\n" );
 312                 fprintf( stderr, "before align 1 %s \n", indication1 );
 313                 display( mseq1, clus1 );
 314                 fprintf( stderr, "\n" );
 315                 fprintf( stderr, "before align 2 %s \n", indication2 );
 316                 display( mseq2, clus2 );
 317                 fprintf( stderr, "\n" );
 318 */
 319
 320                 if( constraint == 2 )
 321                 {
 322                         imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 323                         pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &dumfl );
 324                 }
 325                 else if( use_fft )
 326                 {
 327                         pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 328                 }
 329                 else
 330                 {
 331                         switch( alg )
 332                         {
 333                                 case( 'a' ):
 334                                         pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 335                                         break;
 336                                 case( 'A' ):
 337                                         pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &dumfl );
 338                                         break;
 339                                 case ( 'C' ):
 340                                         if( outgap && ( clus1 == 1 && clus2 != 1 || clus1 != 1 && clus2 == 1 ) )
 341                                         {
 342                                                 pscore = translate_and_Calign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 343                                         }
 344                                         else
 345                                         {
 346                                                 pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &dumfl );
 347                                         }
 348                                         break;
 349                                 default:
 350                                         ErrorExit( "ERROR IN SOURCE FILE" );
 351                         }
 352                 }
 353 #if SCOREOUT
 354                 fprintf( stderr, "score = %10.2f\n", pscore );
 355 #endif
 356                 tscore += pscore;
 357 /*
 358                 fprintf( stderr, "after align 1 %s \n", indication1 );
 359                 display( mseq1, clus1 );
 360                 fprintf( stderr, "\n" );
 361                 fprintf( stderr, "after align 2 %s \n", indication2 );
 362                 display( mseq2, clus2 );
 363                 fprintf( stderr, "\n" );
 364 */
 365
 366                 writePre( njob, name, nlen, aseq, 0 );
 367
 368                 if( disp ) display( aseq, njob );
 369         }
 370 #if SCOREOUT
 371         fprintf( stderr, "totalscore = %10.2f\n\n", tscore );
 372 #endif
 373 }
 374
 375 static void WriteOptions( FILE *fp )
 376 {
 377
 378         if( dorp == 'd' ) fprintf( fp, "DNA\n" );
 379         else
 380         {
 381                 if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
 382                 else if( scoremtx ==  1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
 383                 else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
 384         }
 385     fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 386     if( use_fft ) fprintf( fp, "FFT on\n" );
 387
 388         fprintf( fp, "tree-base method\n" );
 389         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
 390         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
 391         if( tbitr || tbweight )
 392         {
 393                 fprintf( fp, "iterate at each step\n" );
 394                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
 395                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
 396                 if( tbweight ) fprintf( fp, "  weighted\n" );
 397                 fprintf( fp, "\n" );
 398         }
 399
 400          fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 401
 402         if( alg == 'a' )
 403                 fprintf( fp, "Algorithm A\n" );
 404         else if( alg == 'A' )
 405                 fprintf( fp, "Algorithm A+\n" );
 406         else if( alg == 'S' )
 407                 fprintf( fp, "Apgorithm S\n" );
 408         else if( alg == 'C' )
 409                 fprintf( fp, "Apgorithm A+/C\n" );
 410         else
 411                 fprintf( fp, "Unknown algorithm\n" );
 412
 413         if( treemethod == 'x' )
 414                 fprintf( fp, "Tree = UPGMA (3).\n" );
 415         else if( treemethod == 's' )
 416                 fprintf( fp, "Tree = UPGMA (2).\n" );
 417         else if( treemethod == 'p' )
 418                 fprintf( fp, "Tree = UPGMA (1).\n" );
 419         else
 420                 fprintf( fp, "Unknown tree.\n" );
 421
 422     if( use_fft )
 423     {
 424         fprintf( fp, "FFT on\n" );
 425         if( dorp == 'd' )
 426             fprintf( fp, "Basis : 4 nucleotides\n" );
 427         else
 428         {
 429             if( fftscore )
 430                 fprintf( fp, "Basis : Polarity and Volume\n" );
 431             else
 432                 fprintf( fp, "Basis : 20 amino acids\n" );
 433         }
 434         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
 435         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
 436     }
 437         else
 438         fprintf( fp, "FFT off\n" );
 439         fflush( fp );
 440 }
 441
 442
 443 int main( int argc, char *argv[] )
 444 {
 445         static int  nlen[M];
 446         static char name[M][B], **seq;
 447         static char **mseq1, **mseq2;
 448         static char **aseq;
 449         static char **bseq;
 450         static double **iscore;
 451         static double *eff;
 452         int i, j;
 453         int identity;
 454         static int ***topol;
 455         static double **len;
 456         FILE *prep;
 457         FILE *infp;
 458         FILE *hat3p;
 459         FILE *orderfp;
 460         char c;
 461         int alloclen;
 462         double total;
 463         LocalHom **localhomtable;
 464         LocalHom *tmpptr;
 465
 466         arguments( argc, argv );
 467
 468         if( inputfile )
 469         {
 470                 infp = fopen( inputfile, "r" );
 471                 if( !infp )
 472                 {
 473                         fprintf( stderr, "Cannot open %s\n", inputfile );
 474                         exit( 1 );
 475                 }
 476         }
 477         else
 478                 infp = stdin;
 479
 480
 481         getnumlen( infp );
 482         rewind( infp );
 483
 484         if( njob < 2 )
 485         {
 486                 fprintf( stderr, "At least 2 sequences should be input!\n"
 487                                                  "Only %d sequence found.\n", njob );
 488                 exit( 1 );
 489         }
 490
 491         seq = AllocateCharMtx( njob, nlenmax*9+1 );
 492         aseq = AllocateCharMtx( njob, nlenmax*9+1 );
 493         bseq = AllocateCharMtx( njob, nlenmax*9+1 );
 494         mseq1 = AllocateCharMtx( njob, 0 );
 495         mseq2 = AllocateCharMtx( njob, 0 );
 496         alloclen = nlenmax*9;
 497
 498         topol = AllocateIntCub( njob, 2, njob );
 499         len = AllocateDoubleMtx( njob, 2 );
 500         iscore = AllocateDoubleMtx( njob, njob );
 501         eff = AllocateDoubleVec( njob );
 502         if( constraint )
 503         {
 504                 localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 505                 for( i=0; i<njob; i++)
 506                 {
 507                         localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
 508                         for( j=0; j<njob; j++)
 509                         {
 510                                 localhomtable[i][j].nokori = 0;
 511                                 localhomtable[i][j].start1 = -1;
 512                                 localhomtable[i][j].end1 = -1;
 513                                 localhomtable[i][j].start2 = -1;
 514                                 localhomtable[i][j].end2 = -1;
 515                                 localhomtable[i][j].overlapaa = -1.0;
 516                                 localhomtable[i][j].opt = -1.0;
 517                                 localhomtable[i][j].importance = -1.0;
 518                                 localhomtable[i][j].next = NULL;
 519                                 localhomtable[i][j].extended = -1;
 520                                 localhomtable[i][j].last = localhomtable[i]+j;
 521                         }
 522                 }
 523
 524                 fprintf( stderr, "Loading 'hat3' ... " );
 525                 prep = fopen( "hat3", "r" );
 526                 if( prep == NULL ) ErrorExit( "Make hat3." );
 527                 readlocalhomtable2( prep, njob, localhomtable );
 528                 fclose( prep );
 529                 fprintf( stderr, "done.\n" );
 530
 531 #if 0
 532                 for( i=0; i<njob; i++ )
 533                         for( j=0; j<njob; j++ )
 534                                 fprintf( stderr, "[%d][%d].nokori = %d\n", i, j, localhomtable[i][j].nokori );
 535 #endif
 536
 537 //              outlocalhom( localhomtable, njob );
 538
 539         }
 540
 541
 542
 543
 544 #if 0
 545         Read( name, nlen, seq );
 546 #else
 547         readData( infp, name, nlen, seq );
 548 #endif
 549         fclose( infp );
 550
 551         constants( njob, seq );
 552
 553 #if 0
 554         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
 555 #endif
 556
 557         initSignalSM();
 558
 559         initFiles();
 560
 561         WriteOptions( trap_g );
 562
 563         c = seqcheck( seq );
 564         if( c )
 565         {
 566                 fprintf( stderr, "Illeagal character %c\n", c );
 567                 exit( 1 );
 568         }
 569
 570         writePre( njob, name, nlen, seq, 0 );
 571
 572         if( tbutree == 0 )
 573         {
 574                 for( i=1; i<njob; i++ )
 575                 {
 576                         if( nlen[i] != nlen[0] )
 577                         {
 578                                 fprintf( stderr, "Input pre-aligned seqences or make hat2.\n" );
 579                                 exit( 1 );
 580                         }
 581                 }
 582                 for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 583                 {
 584                 /*
 585                         pscore[i][j] = (double)score_calc1( seq[i], seq[j] );
 586                 */
 587                         iscore[i][j] = (double)( substitution_hosei( seq[i], seq[j] ) );
 588                 }
 589         }
 590         else
 591         {
 592                 fprintf( stderr, "Loading 'hat2' ... " );
 593                 prep = fopen( "hat2", "r" );
 594                 if( prep == NULL ) ErrorExit( "Make hat2." );
 595                 readhat2( prep, njob, name, iscore );
 596                 fclose( prep );
 597                 fprintf( stderr, "done.\n" );
 598
 599 #if 0
 600                 prep = fopen( "hat2_check", "w" );
 601                 WriteHat2( prep, njob, name, pscore );
 602                 fclose( prep );
 603 #endif
 604
 605         }
 606
 607 #if 0
 608         for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 609                 iscore[j][i] = iscore[i][j];
 610 #endif
 611
 612         fprintf( stderr, "Constructing dendrogram ... " );
 613         if( treemethod == 'x' )
 614                 veryfastsupg( njob, iscore, topol, len );
 615         else if( treemethod == 's' )
 616                 spg( njob, iscore, topol, len );
 617         else if( treemethod == 'p' )
 618                 upg2( njob, iscore, topol, len );
 619         else
 620                 ErrorExit( "Incorrect tree\n" );
 621         fprintf( stderr, "done.\n" );
 622
 623         orderfp = fopen( "order", "w" );
 624         if( !orderfp )
 625         {
 626                 fprintf( stderr, "Cannot open 'order'\n" );
 627                 exit( 1 );
 628         }
 629         for( i=0; (j=topol[njob-2][0][i])!=-1; i++ )
 630         {
 631                 fprintf( orderfp, "%d\n", j );
 632         }
 633         for( i=0; (j=topol[njob-2][1][i])!=-1; i++ )
 634         {
 635                 fprintf( orderfp, "%d\n", j );
 636         }
 637         fclose( orderfp );
 638
 639
 640 //      countnode( njob, topol, node0 );
 641         if( tbrweight )
 642         {
 643                 weight = 3;
 644 #if 0
 645                 utree = 0; counteff( njob, topol, len, eff ); utree = 1;
 646 #else
 647                 counteff_simple( njob, topol, len, eff );
 648 #endif
 649         }
 650         else
 651         {
 652                 for( i=0; i<njob; i++ ) eff[i] = 1.0;
 653         }
 654
 655         if( constraint )
 656         {
 657 #if 1
 658                 if( thrinter > 0.0 )
 659                 {
 660                         fprintf( stderr, "Extending localhom ... " );
 661                         extendlocalhom2( njob, localhomtable, iscore );
 662                         fprintf( stderr, "done.\n" );
 663 // job de watteoku
 664                 for( i=0; i<njob-1; i++ )
 665                         {
 666                                 for( j=i+1; j<njob; j++ )
 667                                 {
 668                                         if( i == j ) continue;
 669                                         if( localhomtable[i][j]->opt == -1.0 ) continue;
 670                                         for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
 671                                                 tmpptr->opt /= njob;
 672                                 }
 673                         }
 674                 }
 675 #endif
 676 //              resetlocalhom( njob, localhomtable ); // chuui !!!
 677 #if 1
 678             fprintf( stderr, "##### writing hat4, nlenmax = %d\n", nlenmax );
 679         hat3p = fopen( "hat4", "w" );
 680         if( !hat3p ) ErrorExit( "Cannot open hat4." );
 681         for( i=0; i<njob-1; i++ )
 682         {
 683                 for( j=i+1; j<njob; j++ )
 684                 {
 685                                 if( i == j ) continue;
 686                 for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
 687                 {
 688                         fprintf( hat3p, "%d %d %d %6.3f %d %d %d %d %p nokori=%d ext=%d\n", i, j, tmpptr->overlapaa, tmpptr->opt * 5.8 / 600, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->next, tmpptr->nokori, tmpptr->extended );
 689                                         fflush( hat3p );
 690                         if( tmpptr->opt == -1.0 ) continue;
 691                 }
 692                 }
 693         }
 694         fclose( hat3p );
 695             fprintf( stderr, "done.\n" );
 696 #endif
 697         }
 698
 699         FreeDoubleMtx( len );
 700         FreeDoubleMtx( iscore );
 701
 702         for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
 703
 704         treebase( name, nlen, bseq, aseq, mseq1, mseq2, topol, eff, alloclen, localhomtable );
 705
 706
 707 #if 0
 708         if( constraint )
 709         {
 710                 LocalHom *tmppt1, *tmppt2;
 711                 for( i=0; i<njob; i++)
 712                 {
 713                         for( j=0; j<njob; j++)
 714                         {
 715                                 tmppt1 = localhomtable[i]+j;
 716                                 while( tmppt2 = tmppt1->next )
 717                                 {
 718                                         free( (void *)tmppt1 );
 719                                         tmppt1 = tmppt2;
 720                                 }
 721                                 free( (void *)tmppt1 );
 722                         }
 723                         free( (void *)(localhomtable[i]+j) );
 724                 }
 725                 free( (void *)localhomtable );
 726         }
 727 #endif
 728
 729         fprintf( trap_g, "done.\n" );
 730         fclose( trap_g );
 731
 732         writePre( njob, name, nlen, aseq, !contin );
 733 #if 0
 734         writeData( stdout, njob, name, nlen, aseq );
 735 #endif
 736 #if IODEBUG
 737         fprintf( stderr, "OSHIMAI\n" );
 738 #endif
 739
 740         if( constraint ) FreeLocalHomTable( localhomtable, njob );
 741
 742         SHOWVERSION;
 743         return( 0 );
 744 }