binaries/src/mafft/core/pairlocalalign.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 0
   6
   7 #define NODIST -9999
   8
   9 static char *whereispairalign;
  10 static char *laraparams;
  11 static char foldalignopt[1000];
  12 static int stdout_align;
  13 static int stdout_dist;
  14 static int store_localhom;
  15 static int store_dist;
  16
  17 typedef struct _jobtable
  18 {
  19         int i;
  20         int j;
  21 } Jobtable;
  22
  23 #ifdef enablemultithread
  24 typedef struct _thread_arg
  25 {
  26         int thread_no;
  27         int njob;
  28         Jobtable *jobpospt;
  29         char **name;
  30         char **seq;
  31         LocalHom **localhomtable;
  32         double **distancemtx;
  33         double *selfscore;
  34         char ***bpp;
  35         int alloclen;
  36         pthread_mutex_t *mutex_counter;
  37         pthread_mutex_t *mutex_stdout;
  38 } thread_arg_t;
  39 #endif
  40
  41 static void t2u( char *seq )
  42 {
  43         while( *seq )
  44         {
  45                 if     ( *seq == 'A' ) *seq = 'a';
  46                 else if( *seq == 'a' ) *seq = 'a';
  47                 else if( *seq == 'T' ) *seq = 'u';
  48                 else if( *seq == 't' ) *seq = 'u';
  49                 else if( *seq == 'U' ) *seq = 'u';
  50                 else if( *seq == 'u' ) *seq = 'u';
  51                 else if( *seq == 'G' ) *seq = 'g';
  52                 else if( *seq == 'g' ) *seq = 'g';
  53                 else if( *seq == 'C' ) *seq = 'c';
  54                 else if( *seq == 'c' ) *seq = 'c';
  55                 else *seq = 'n';
  56                 seq++;
  57         }
  58 }
  59
  60 static float recallpairfoldalign( char **mseq1, char **mseq2, int m1, int m2, int *of1pt, int *of2pt, int alloclen )
  61 {
  62         static FILE *fp = NULL;
  63         float value;
  64         char *aln1;
  65         char *aln2;
  66         int of1tmp, of2tmp;
  67
  68         if( fp == NULL )
  69         {
  70                 fp = fopen( "_foldalignout", "r" );
  71                 if( fp == NULL )
  72                 {
  73                         fprintf( stderr, "Cannot open _foldalignout\n" );
  74                         exit( 1 );
  75                 }
  76         }
  77
  78         aln1 = calloc( alloclen, sizeof( char ) );
  79         aln2 = calloc( alloclen, sizeof( char ) );
  80
  81         readpairfoldalign( fp, *mseq1, *mseq2, aln1, aln2, m1, m2, &of1tmp, &of2tmp, alloclen );
  82
  83         if( strstr( foldalignopt, "-global") )
  84         {
  85                 fprintf( stderr, "Calling G__align11\n" );
  86                 value = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
  87                 *of1pt = 0;
  88                 *of2pt = 0;
  89         }
  90         else
  91         {
  92                 fprintf( stderr, "Calling L__align11\n" );
  93                 value = L__align11( mseq1, mseq2, alloclen, of1pt, of2pt );
  94         }
  95
  96 //      value = (float)naivepairscore11( *mseq1, *mseq2, penalty ); // nennnotame
  97
  98         if( aln1[0] == 0 )
  99         {
 100                 fprintf( stderr, "FOLDALIGN returned no alignment between %d and %d.  Sequence alignment is used instead.\n", m1+1, m2+1 );
 101         }
 102         else
 103         {
 104                 strcpy( *mseq1, aln1 );
 105                 strcpy( *mseq2, aln2 );
 106                 *of1pt = of1tmp;
 107                 *of2pt = of2tmp;
 108         }
 109
 110 //      value = naivepairscore11( *mseq1, *mseq2, penalty ); // v6.511 ha kore wo tsukau, global nomi dakara.
 111
 112 //      fclose( fp ); // saigo dake yatta houga yoi.
 113
 114 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 115 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 116
 117
 118         free( aln1 );
 119         free( aln2 );
 120
 121         return( value );
 122 }
 123
 124 static void callfoldalign( int nseq, char **mseq )
 125 {
 126         FILE *fp;
 127         int i;
 128         int res;
 129         static char com[10000];
 130
 131         for( i=0; i<nseq; i++ )
 132                 t2u( mseq[i] );
 133
 134         fp = fopen( "_foldalignin", "w" );
 135         if( !fp )
 136         {
 137                 fprintf( stderr, "Cannot open _foldalignin\n" );
 138                 exit( 1 );
 139         }
 140         for( i=0; i<nseq; i++ )
 141         {
 142                 fprintf( fp, ">%d\n", i+1 );
 143                 fprintf( fp, "%s\n", mseq[i] );
 144         }
 145         fclose( fp );
 146
 147         sprintf( com, "env PATH=%s  foldalign210 %s _foldalignin > _foldalignout ", whereispairalign, foldalignopt );
 148         res = system( com );
 149         if( res )
 150         {
 151                 fprintf( stderr, "Error in foldalign\n" );
 152                 exit( 1 );
 153         }
 154
 155 }
 156
 157 static void calllara( int nseq, char **mseq, char *laraarg )
 158 {
 159         FILE *fp;
 160         int i;
 161         int res;
 162         static char com[10000];
 163
 164 //      for( i=0; i<nseq; i++ )
 165
 166         fp = fopen( "_larain", "w" );
 167         if( !fp )
 168         {
 169                 fprintf( stderr, "Cannot open _larain\n" );
 170                 exit( 1 );
 171         }
 172         for( i=0; i<nseq; i++ )
 173         {
 174                 fprintf( fp, ">%d\n", i+1 );
 175                 fprintf( fp, "%s\n", mseq[i] );
 176         }
 177         fclose( fp );
 178
 179
 180 //      fprintf( stderr, "calling LaRA\n" );
 181         sprintf( com, "env PATH=%s:/bin:/usr/bin mafft_lara -i _larain -w _laraout -o _lara.params %s", whereispairalign, laraarg );
 182         res = system( com );
 183         if( res )
 184         {
 185                 fprintf( stderr, "Error in lara\n" );
 186                 exit( 1 );
 187         }
 188 }
 189
 190 static float recalllara( char **mseq1, char **mseq2, int alloclen )
 191 {
 192         static FILE *fp = NULL;
 193         static char *ungap1;
 194         static char *ungap2;
 195         static char *ori1;
 196         static char *ori2;
 197 //      int res;
 198         static char com[10000];
 199         float value;
 200
 201
 202         if( fp == NULL )
 203         {
 204                 fp = fopen( "_laraout", "r" );
 205                 if( fp == NULL )
 206                 {
 207                         fprintf( stderr, "Cannot open _laraout\n" );
 208                         exit( 1 );
 209                 }
 210                 ungap1 = AllocateCharVec( alloclen );
 211                 ungap2 = AllocateCharVec( alloclen );
 212                 ori1 = AllocateCharVec( alloclen );
 213                 ori2 = AllocateCharVec( alloclen );
 214         }
 215
 216
 217         strcpy( ori1, *mseq1 );
 218         strcpy( ori2, *mseq2 );
 219
 220         fgets( com, 999, fp );
 221         myfgets( com, 9999, fp );
 222         strcpy( *mseq1, com );
 223         myfgets( com, 9999, fp );
 224         strcpy( *mseq2, com );
 225
 226         gappick0( ungap1, *mseq1 );
 227         gappick0( ungap2, *mseq2 );
 228         t2u( ungap1 );
 229         t2u( ungap2 );
 230         t2u( ori1 );
 231         t2u( ori2 );
 232
 233         if( strcmp( ungap1, ori1 ) || strcmp( ungap2, ori2 ) )
 234         {
 235                 fprintf( stderr, "SEQUENCE CHANGED!!\n" );
 236                 fprintf( stderr, "*mseq1  = %s\n", *mseq1 );
 237                 fprintf( stderr, "ungap1  = %s\n", ungap1 );
 238                 fprintf( stderr, "ori1    = %s\n", ori1 );
 239                 fprintf( stderr, "*mseq2  = %s\n", *mseq2 );
 240                 fprintf( stderr, "ungap2  = %s\n", ungap2 );
 241                 fprintf( stderr, "ori2    = %s\n", ori2 );
 242                 exit( 1 );
 243         }
 244
 245         value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
 246
 247 //      fclose( fp ); // saigo dake yatta houga yoi.
 248
 249         return( value );
 250 }
 251
 252
 253 static float callmxscarna_giving_bpp( char **mseq1, char **mseq2, char **bpp1, char **bpp2, int alloclen, int i, int j )
 254 {
 255         FILE *fp;
 256         int res;
 257         char *com;
 258         float value;
 259         char *dirname;
 260
 261
 262         dirname = calloc( 100, sizeof( char ) );
 263         com = calloc( 1000, sizeof( char ) );
 264         sprintf( dirname, "_%d-%d", i, j );
 265         sprintf( com, "rm -rf %s", dirname );
 266         system( com );
 267         sprintf( com, "mkdir %s", dirname );
 268         system( com );
 269
 270
 271         sprintf( com, "%s/_bpporg", dirname );
 272         fp = fopen( com, "w" );
 273         if( !fp )
 274         {
 275                 fprintf( stderr, "Cannot write to %s/_bpporg\n", dirname );
 276                 exit( 1 );
 277         }
 278         fprintf( fp, ">a\n" );
 279         while( *bpp1 )
 280                 fprintf( fp, "%s", *bpp1++ );
 281
 282         fprintf( fp, ">b\n" );
 283         while( *bpp2 )
 284                 fprintf( fp, "%s", *bpp2++ );
 285         fclose( fp );
 286
 287         sprintf( com, "tr -d '\\r' < %s/_bpporg > %s/_bpp", dirname, dirname );
 288         system( com ); // for cygwin, wakaran
 289
 290         t2u( *mseq1 );
 291         t2u( *mseq2 );
 292
 293         sprintf( com, "%s/_mxscarnainorg", dirname );
 294         fp = fopen( com, "w" );
 295         if( !fp )
 296         {
 297                 fprintf( stderr, "Cannot open %s/_mxscarnainorg\n", dirname );
 298                 exit( 1 );
 299         }
 300         fprintf( fp, ">1\n" );
 301 //      fprintf( fp, "%s\n", *mseq1 );
 302         write1seq( fp, *mseq1 );
 303         fprintf( fp, ">2\n" );
 304 //      fprintf( fp, "%s\n", *mseq2 );
 305         write1seq( fp, *mseq2 );
 306         fclose( fp );
 307
 308         sprintf( com, "tr -d '\\r' < %s/_mxscarnainorg > %s/_mxscarnain", dirname, dirname );
 309         system( com ); // for cygwin, wakaran
 310
 311 #if 0
 312         sprintf( com, "cd %s; %s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum", dirname, whereispairalign );
 313 #else
 314         sprintf( com, "_mxscarnash%s", dirname );
 315         fp = fopen( com, "w" );
 316         fprintf( fp, "cd %s\n", dirname );
 317         fprintf( fp, "%s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum\n", whereispairalign );
 318         fprintf( fp, "exit $tatus\n" );
 319         fclose( fp );
 320
 321         sprintf( com, "tr -d '\\r' < _mxscarnash%s > _mxscarnash%s.unix", dirname, dirname );
 322         system( com ); // for cygwin, wakaran
 323
 324         sprintf( com, "sh _mxscarnash%s.unix 2>_dum%s", dirname, dirname );
 325 #endif
 326         res = system( com );
 327         if( res )
 328         {
 329                 fprintf( stderr, "Error in mxscarna\n" );
 330                 exit( 1 );
 331         }
 332
 333         sprintf( com, "%s/_mxscarnaout", dirname );
 334
 335         fp = fopen( com, "r" );
 336         if( !fp )
 337         {
 338                 fprintf( stderr, "Cannot open %s/_mxscarnaout\n", dirname );
 339                 exit( 1 );
 340         }
 341
 342         fgets( com, 999, fp );
 343         load1SeqWithoutName_new( fp, *mseq1 );
 344         fgets( com, 999, fp );
 345         load1SeqWithoutName_new( fp, *mseq2 );
 346
 347         fclose( fp );
 348
 349 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 350 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 351
 352         value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
 353
 354 #if 0
 355         sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
 356         if( system( com ) )
 357         {
 358                 fprintf( stderr, "retrying to rmdir\n" );
 359                 usleep( 2000 );
 360                 system( com );
 361         }
 362 #endif
 363
 364         free( dirname );
 365         free( com );
 366
 367
 368         return( value );
 369 }
 370
 371 #if 0
 372 static float callmxscarna_slow( char **mseq1, char **mseq2, int alloclen )
 373 {
 374         FILE *fp;
 375         int res;
 376         static char com[10000];
 377         float value;
 378
 379
 380         t2u( *mseq1 );
 381         t2u( *mseq2 );
 382         fp = fopen( "_mxscarnain", "w" );
 383         if( !fp )
 384         {
 385                 fprintf( stderr, "Cannot open _mxscarnain\n" );
 386                 exit( 1 );
 387         }
 388         fprintf( fp, ">1\n" );
 389         fprintf( fp, "%s\n", *mseq1 );
 390         fprintf( fp, ">2\n" );
 391         fprintf( fp, "%s\n", *mseq2 );
 392         fclose( fp );
 393
 394         sprintf( com, "env PATH=%s mxscarnamod _mxscarnain > _mxscarnaout 2>_dum", whereispairalign );
 395         res = system( com );
 396         if( res )
 397         {
 398                 fprintf( stderr, "Error in mxscarna\n" );
 399                 exit( 1 );
 400         }
 401
 402         fp = fopen( "_mxscarnaout", "r" );
 403         if( !fp )
 404         {
 405                 fprintf( stderr, "Cannot open _mxscarnaout\n" );
 406                 exit( 1 );
 407         }
 408
 409         fgets( com, 999, fp );
 410         load1SeqWithoutName_new( fp, *mseq1 );
 411         fgets( com, 999, fp );
 412         load1SeqWithoutName_new( fp, *mseq2 );
 413
 414         fclose( fp );
 415
 416 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 417 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 418
 419         value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
 420
 421         return( value );
 422 }
 423 #endif
 424
 425 static void readhat4( FILE *fp, char ***bpp )
 426 {
 427         char oneline[1000];
 428         int bppsize;
 429         int onechar;
 430 //      double prob;
 431 //      int posi, posj;
 432
 433         bppsize = 0;
 434 //      fprintf( stderr, "reading hat4\n" );
 435         onechar = getc(fp);
 436 //      fprintf( stderr, "onechar = %c\n", onechar );
 437         if( onechar != '>' )
 438         {
 439                 fprintf( stderr, "Format error\n" );
 440                 exit( 1 );
 441         }
 442         ungetc( onechar, fp );
 443         fgets( oneline, 999, fp );
 444         while( 1 )
 445         {
 446                 onechar = getc(fp);
 447                 ungetc( onechar, fp );
 448                 if( onechar == '>' || onechar == EOF )
 449                 {
 450 //                      fprintf( stderr, "Next\n" );
 451                         *bpp = realloc( *bpp, (bppsize+2) * sizeof( char * ) );
 452                         (*bpp)[bppsize] = NULL;
 453                         break;
 454                 }
 455                 fgets( oneline, 999, fp );
 456 //              fprintf( stderr, "oneline=%s\n", oneline );
 457 //              sscanf( oneline, "%d %d %f", &posi, &posj, &prob );
 458 //              fprintf( stderr, "%d %d -> %f\n", posi, posj, prob );
 459                 *bpp = realloc( *bpp, (bppsize+2) * sizeof( char * ) );
 460                 (*bpp)[bppsize] = calloc( 100, sizeof( char ) );
 461                 strcpy( (*bpp)[bppsize], oneline );
 462                 bppsize++;
 463         }
 464 }
 465
 466 static void preparebpp( int nseq, char ***bpp )
 467 {
 468         FILE *fp;
 469         int i;
 470
 471         fp = fopen( "hat4", "r" );
 472         if( !fp )
 473         {
 474                 fprintf( stderr, "Cannot open hat4\n" );
 475                 exit( 1 );
 476         }
 477         for( i=0; i<nseq; i++ )
 478                 readhat4( fp, bpp+i );
 479         fclose( fp );
 480 }
 481
 482 void arguments( int argc, char *argv[] )
 483 {
 484     int c;
 485
 486         nthread = 1;
 487         foldalignopt[0] = 0;
 488         laraparams = NULL;
 489         inputfile = NULL;
 490         fftkeika = 0;
 491         pslocal = -1000.0;
 492         constraint = 0;
 493         nblosum = 62;
 494         fmodel = 0;
 495         calledByXced = 0;
 496         devide = 0;
 497         use_fft = 0;
 498         fftscore = 1;
 499         fftRepeatStop = 0;
 500         fftNoAnchStop = 0;
 501     weight = 3;
 502     utree = 1;
 503         tbutree = 1;
 504     refine = 0;
 505     check = 1;
 506     cut = 0.0;
 507     disp = 0;
 508     outgap = 1;
 509     alg = 'A';
 510     mix = 0;
 511         tbitr = 0;
 512         scmtd = 5;
 513         tbweight = 0;
 514         tbrweight = 3;
 515         checkC = 0;
 516         treemethod = 'x';
 517         contin = 0;
 518         scoremtx = 1;
 519         kobetsubunkatsu = 0;
 520         divpairscore = 0;
 521         stdout_align = 0;
 522         stdout_dist = 0;
 523         store_dist = 1;
 524         store_localhom = 1;
 525         dorp = NOTSPECIFIED;
 526         ppenalty = NOTSPECIFIED;
 527         ppenalty_OP = NOTSPECIFIED;
 528         ppenalty_ex = NOTSPECIFIED;
 529         ppenalty_EX = NOTSPECIFIED;
 530         poffset = NOTSPECIFIED;
 531         kimuraR = NOTSPECIFIED;
 532         pamN = NOTSPECIFIED;
 533         geta2 = GETA2;
 534         fftWinSize = NOTSPECIFIED;
 535         fftThreshold = NOTSPECIFIED;
 536         RNAppenalty = NOTSPECIFIED;
 537         RNApthr = NOTSPECIFIED;
 538
 539     while( --argc > 0 && (*++argv)[0] == '-' )
 540         {
 541         while ( ( c = *++argv[0] ) )
 542                 {
 543             switch( c )
 544             {
 545                                 case 'i':
 546                                         inputfile = *++argv;
 547                                         fprintf( stderr, "inputfile = %s\n", inputfile );
 548                                         --argc;
 549                                         goto nextoption;
 550                                 case 'f':
 551                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 552                                         --argc;
 553                                         goto nextoption;
 554                                 case 'g':
 555                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
 556                                         --argc;
 557                                         goto nextoption;
 558                                 case 'O':
 559                                         ppenalty_OP = (int)( atof( *++argv ) * 1000 - 0.5 );
 560                                         --argc;
 561                                         goto nextoption;
 562                                 case 'E':
 563                                         ppenalty_EX = (int)( atof( *++argv ) * 1000 - 0.5 );
 564                                         --argc;
 565                                         goto nextoption;
 566                                 case 'h':
 567                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
 568                                         --argc;
 569                                         goto nextoption;
 570                                 case 'k':
 571                                         kimuraR = atoi( *++argv );
 572 //                                      fprintf( stderr, "kimuraR = %d\n", kimuraR );
 573                                         --argc;
 574                                         goto nextoption;
 575                                 case 'b':
 576                                         nblosum = atoi( *++argv );
 577                                         scoremtx = 1;
 578 //                                      fprintf( stderr, "blosum %d\n", nblosum );
 579                                         --argc;
 580                                         goto nextoption;
 581                                 case 'j':
 582                                         pamN = atoi( *++argv );
 583                                         scoremtx = 0;
 584                                         TMorJTT = JTT;
 585                                         fprintf( stderr, "jtt %d\n", pamN );
 586                                         --argc;
 587                                         goto nextoption;
 588                                 case 'm':
 589                                         pamN = atoi( *++argv );
 590                                         scoremtx = 0;
 591                                         TMorJTT = TM;
 592                                         fprintf( stderr, "TM %d\n", pamN );
 593                                         --argc;
 594                                         goto nextoption;
 595                                 case 'l':
 596                                         ppslocal = (int)( atof( *++argv ) * 1000 + 0.5 );
 597                                         pslocal = (int)( 600.0 / 1000.0 * ppslocal + 0.5);
 598 //                                      fprintf( stderr, "ppslocal = %d\n", ppslocal );
 599 //                                      fprintf( stderr, "pslocal = %d\n", pslocal );
 600                                         --argc;
 601                                         goto nextoption;
 602                                 case 'd':
 603                                         whereispairalign = *++argv;
 604                                         fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
 605                                         --argc;
 606                                         goto nextoption;
 607                                 case 'p':
 608                                         laraparams = *++argv;
 609                                         fprintf( stderr, "laraparams = %s\n", laraparams );
 610                                         --argc;
 611                                         goto nextoption;
 612                                 case 'C':
 613                                         nthread = atoi( *++argv );
 614                                         fprintf( stderr, "nthread = %d\n", nthread );
 615                                         --argc;
 616                                         goto nextoption;
 617                                 case 'c':
 618                                         stdout_dist = 1;
 619                                         break;
 620                                 case 'n':
 621                                         stdout_align = 1;
 622                                         break;
 623                                 case 'x':
 624                                         store_localhom = 0;
 625                                         store_dist = 0;
 626                                         break;
 627 #if 1
 628                                 case 'a':
 629                                         fmodel = 1;
 630                                         break;
 631 #endif
 632                                 case 'r':
 633                                         fmodel = -1;
 634                                         break;
 635                                 case 'D':
 636                                         dorp = 'd';
 637                                         break;
 638                                 case 'P':
 639                                         dorp = 'p';
 640                                         break;
 641                                 case 'e':
 642                                         fftscore = 0;
 643                                         break;
 644 #if 0
 645                                 case 'O':
 646                                         fftNoAnchStop = 1;
 647                                         break;
 648 #endif
 649                                 case 'Q':
 650                                         calledByXced = 1;
 651                                         break;
 652 #if 0
 653                                 case 'x':
 654                                         disp = 1;
 655                                         break;
 656                                 case 'a':
 657                                         alg = 'a';
 658                                         break;
 659 #endif
 660                                 case 'S':
 661                                         alg = 'S';
 662                                         break;
 663                                 case 't':
 664                                         alg = 't';
 665                                         store_localhom = 0;
 666                                         break;
 667                                 case 'L':
 668                                         alg = 'L';
 669                                         break;
 670                                 case 's':
 671                                         alg = 's';
 672                                         break;
 673                                 case 'B':
 674                                         alg = 'B';
 675                                         break;
 676                                 case 'T':
 677                                         alg = 'T';
 678                                         break;
 679                                 case 'H':
 680                                         alg = 'H';
 681                                         break;
 682                                 case 'M':
 683                                         alg = 'M';
 684                                         break;
 685                                 case 'R':
 686                                         alg = 'R';
 687                                         break;
 688                                 case 'N':
 689                                         alg = 'N';
 690                                         break;
 691                                 case 'K':
 692                                         alg = 'K';
 693                                         break;
 694                                 case 'A':
 695                                         alg = 'A';
 696                                         break;
 697                                 case 'V':
 698                                         alg = 'V';
 699                                         break;
 700                                 case 'F':
 701                                         use_fft = 1;
 702                                         break;
 703                                 case 'v':
 704                                         tbrweight = 3;
 705                                         break;
 706                                 case 'y':
 707                                         divpairscore = 1;
 708                                         break;
 709 /* Modified 01/08/27, default: user tree */
 710                                 case 'J':
 711                                         tbutree = 0;
 712                                         break;
 713 /* modification end. */
 714                                 case 'o':
 715 //                                      foldalignopt = *++argv;
 716                                         strcat( foldalignopt, " " );
 717                                         strcat( foldalignopt, *++argv );
 718                                         fprintf( stderr, "foldalignopt = %s\n", foldalignopt );
 719                                         --argc;
 720                                         goto nextoption;
 721                                 case 'z':
 722                                         fftThreshold = atoi( *++argv );
 723                                         --argc;
 724                                         goto nextoption;
 725                                 case 'w':
 726                                         fftWinSize = atoi( *++argv );
 727                                         --argc;
 728                                         goto nextoption;
 729                                 case 'Z':
 730                                         checkC = 1;
 731                                         break;
 732                 default:
 733                     fprintf( stderr, "illegal option %c\n", c );
 734                     argc = 0;
 735                     break;
 736             }
 737                 }
 738                 nextoption:
 739                         ;
 740         }
 741     if( argc == 1 )
 742     {
 743         cut = atof( (*argv) );
 744         argc--;
 745     }
 746     if( argc != 0 )
 747     {
 748         fprintf( stderr, "options: Check source file !\n" );
 749         exit( 1 );
 750     }
 751         if( tbitr == 1 && outgap == 0 )
 752         {
 753                 fprintf( stderr, "conflicting options : o, m or u\n" );
 754                 exit( 1 );
 755         }
 756 }
 757
 758 int countamino( char *s, int end )
 759 {
 760         int val = 0;
 761         while( end-- )
 762                 if( *s++ != '-' ) val++;
 763         return( val );
 764 }
 765
 766 #if enablemultithread
 767 static void *athread( void *arg )
 768 {
 769         thread_arg_t *targ = (thread_arg_t *)arg;
 770         int i, j;
 771         int clus1, clus2;
 772         int off1, off2;
 773         int intdum;
 774         double bunbo;
 775         float pscore = 0.0; // by D.Mathog
 776         double *effarr;
 777         double *effarr1;
 778         double *effarr2;
 779         char **pair;
 780         char *indication1, *indication2;
 781         char **mseq1, **mseq2;
 782         char **aseq;
 783
 784 // thread_arg
 785         int thread_no = targ->thread_no;
 786         int njob = targ->njob;
 787         Jobtable *jobpospt = targ->jobpospt;
 788         char **name = targ->name;
 789         char **seq = targ->seq;
 790         LocalHom **localhomtable = targ->localhomtable;
 791         double **distancemtx = targ->distancemtx;
 792         double *selfscore = targ->selfscore;
 793         char ***bpp = targ->bpp;
 794         int alloclen = targ->alloclen;
 795
 796 //      fprintf( stderr, "thread %d start!\n", thread_no );
 797
 798         effarr = AllocateDoubleVec( njob );
 799         for( i=0; i<njob; i++ ) effarr[i] = 1.0;
 800         effarr1 = AllocateDoubleVec( njob );
 801         effarr2 = AllocateDoubleVec( njob );
 802         indication1 = AllocateCharVec( 150 );
 803         indication2 = AllocateCharVec( 150 );
 804         pair = AllocateCharMtx( njob, njob );
 805         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
 806         for( i=0; i<njob; i++ ) pair[i][i] = 1;
 807         mseq1 = AllocateCharMtx( njob, 0 );
 808         mseq2 = AllocateCharMtx( njob, 0 );
 809         aseq = AllocateCharMtx( njob, alloclen+10 );
 810
 811         while( 1 )
 812         {
 813                 pthread_mutex_lock( targ->mutex_counter );
 814                 j = jobpospt->j;
 815                 i = jobpospt->i;
 816                 j++;
 817                 if( j == njob )
 818                 {
 819                         i++;
 820                         j = i + 1;
 821                         if( i == njob-1 )
 822                         {
 823 //                              fprintf( stderr, "thread %d end!\n", thread_no );
 824                                 pthread_mutex_unlock( targ->mutex_counter );
 825
 826                                 if( commonIP ) FreeIntMtx( commonIP );
 827                                 commonIP = NULL;
 828                                 if( commonJP ) FreeIntMtx( commonJP );
 829                                 commonJP = NULL;
 830                                 Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
 831                                 G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
 832                                 L__align11( NULL, NULL, 0, NULL, NULL );
 833                                 genL__align11( NULL, NULL, 0, NULL, NULL );
 834                                 free( effarr );
 835                                 free( effarr1 );
 836                                 free( effarr2 );
 837                                 free( indication1 );
 838                                 free( indication2 );
 839                                 FreeCharMtx( pair  );
 840                                 free( mseq1 );
 841                                 free( mseq2 );
 842                                 FreeCharMtx( aseq  );
 843                                 return( NULL );
 844                         }
 845                 }
 846                 jobpospt->j = j;
 847                 jobpospt->i = i;
 848                 pthread_mutex_unlock( targ->mutex_counter );
 849
 850
 851                 if( j == i+1 || j % 100 == 0 )
 852                 {
 853                         fprintf( stderr, "% 5d / %d (by thread %3d) \r", i, njob, thread_no );
 854 //                      fprintf( stderr, "% 5d - %5d / %d (thread %d)\n", i, j, njob, thread_no );
 855                 }
 856
 857
 858                 if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
 859                 {
 860                         if( store_dist ) distancemtx[i][j] = 2.0;
 861                         if( stdout_dist)
 862                         {
 863                                 pthread_mutex_lock( targ->mutex_stdout );
 864                                 fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 2.0 );
 865                                 pthread_mutex_unlock( targ->mutex_stdout );
 866                         }
 867                         continue;
 868                 }
 869
 870                 strcpy( aseq[i], seq[i] );
 871                 strcpy( aseq[j], seq[j] );
 872                 clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
 873                 clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
 874 //              fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
 875 //              fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
 876
 877 #if 0
 878                 fprintf( stderr, "group1 = %.66s", indication1 );
 879                 fprintf( stderr, "\n" );
 880                 fprintf( stderr, "group2 = %.66s", indication2 );
 881                 fprintf( stderr, "\n" );
 882 #endif
 883 //              for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
 884
 885                 if( use_fft )
 886                 {
 887                         pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
 888 //                      fprintf( stderr, "pscore (fft) = %f\n", pscore );
 889                         off1 = off2 = 0;
 890                 }
 891                 else
 892                 {
 893                         switch( alg )
 894                         {
 895                                 case( 'L' ):
 896                                         pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
 897                                         L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
 898                                         break;
 899                                 case( 'A' ):
 900                                         pscore = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
 901                                         off1 = off2 = 0;
 902                                         break;
 903                                 case( 'N' ):
 904                                         pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
 905                                         genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
 906                                         break;
 907                                 case( 't' ):
 908                                         pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
 909                                         off1 = off2 = 0;
 910                                         break;
 911                                 case( 's' ):
 912                                         pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
 913                                         off1 = off2 = 0;
 914                                         break;
 915 #if 0
 916                                 case( 'a' ):
 917                                         pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 918                                         off1 = off2 = 0;
 919                                         break;
 920                                 case( 'K' ):
 921                                         pscore = genG__align11( mseq1, mseq2, alloclen );
 922                                         off1 = off2 = 0;
 923                                         break;
 924                                 case( 'H' ):
 925                                         pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
 926                                         break;
 927                                 case( 'B' ):
 928                                 case( 'T' ):
 929                                         pscore = recalllara( mseq1, mseq2, alloclen );
 930                                         off1 = off2 = 0;
 931                                         break;
 932                                 case( 'M' ):
 933                                         pscore = MSalign11( mseq1, mseq2, alloclen );
 934                                         break;
 935 #endif
 936                                 default:
 937                                         ErrorExit( "\n\nERROR IN SOURCE FILE\n\n" );
 938                         }
 939                 }
 940
 941                 if( alg == 't' || ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) ) // 't' no jouken ha iranai to omou. if( ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) )
 942                 {
 943 #if SCOREOUT
 944                         fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
 945 #endif
 946                         if( !store_localhom )
 947                                 ;
 948                         else if( alg == 'H' )
 949                                 putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
 950                         else if( alg != 'S' && alg != 'V' )
 951                         {
 952                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
 953                         }
 954
 955                         if( (bunbo=MIN( selfscore[i], selfscore[j] )) == 0.0 || bunbo < pscore )
 956                                 pscore = 2.0;
 957                         else
 958                                 pscore = ( 1.0 - pscore / bunbo ) * 2.0;
 959                 }
 960                 else
 961                 {
 962                         pscore = 2.0;
 963                 }
 964
 965 #if 1 // mutex
 966                 if( stdout_align )
 967                 {
 968                         pthread_mutex_lock( targ->mutex_stdout );
 969                         if( alg != 't' )
 970                         {
 971                                 fprintf( stdout, "sequence %d - sequence %d, pairwise distance = %10.5f\n", i+1, j+1, pscore );
 972                                 fprintf( stdout, ">%s\n", name[i] );
 973                                 write1seq( stdout, mseq1[0] );
 974                                 fprintf( stdout, ">%s\n", name[j] );
 975                                 write1seq( stdout, mseq2[0] );
 976                                 fprintf( stdout, "\n" );
 977                         }
 978                         pthread_mutex_unlock( targ->mutex_stdout );
 979                 }
 980                 if( stdout_dist )
 981                 {
 982                         pthread_mutex_lock( targ->mutex_stdout );
 983                         if( j == i+1 ) fprintf( stdout, "%d %d d=%.3f\n", i+1, i+1, 0.0 );
 984                         fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, pscore );
 985                         pthread_mutex_unlock( targ->mutex_stdout );
 986                 }
 987 #endif // mutex
 988                 if( store_dist) distancemtx[i][j] = pscore;
 989         }
 990 }
 991 #endif
 992
 993 static void pairalign( char **name, int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, double *effarr, int alloclen )
 994 {
 995         int i, j, ilim;
 996         int clus1, clus2;
 997         int off1, off2;
 998         float pscore = 0.0; // by D.Mathog
 999         static char *indication1, *indication2;
1000         FILE *hat2p, *hat3p;
1001         double **distancemtx;
1002         double *selfscore;
1003         double *effarr1;
1004         double *effarr2;
1005         char *pt;
1006         char *hat2file = "hat2";
1007         LocalHom **localhomtable = NULL, *tmpptr;
1008         static char **pair;
1009         int intdum;
1010         double bunbo;
1011         char ***bpp = NULL; // mxscarna no toki dake
1012
1013         if( store_localhom )
1014         {
1015                 localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
1016                 for( i=0; i<njob; i++)
1017                 {
1018                         localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1019                         for( j=0; j<njob; j++)
1020                         {
1021                                 localhomtable[i][j].start1 = -1;
1022                                 localhomtable[i][j].end1 = -1;
1023                                 localhomtable[i][j].start2 = -1;
1024                                 localhomtable[i][j].end2 = -1;
1025                                 localhomtable[i][j].opt = -1.0;
1026                                 localhomtable[i][j].next = NULL;
1027                                 localhomtable[i][j].nokori = 0;
1028                         }
1029                 }
1030         }
1031
1032         if( store_dist ) distancemtx = AllocateDoubleMtx( njob, njob );
1033         else distancemtx = NULL;
1034         selfscore = AllocateDoubleVec( njob );
1035         effarr1 = AllocateDoubleVec( njob );
1036         effarr2 = AllocateDoubleVec( njob );
1037         indication1 = AllocateCharVec( 150 );
1038         indication2 = AllocateCharVec( 150 );
1039 #if 0
1040 #else
1041         pair = AllocateCharMtx( njob, njob );
1042 #endif
1043
1044 #if 0
1045         fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
1046 #endif
1047
1048 #if 0
1049         for( i=0; i<njob; i++ )
1050                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
1051 #endif
1052
1053
1054 //      writePre( njob, name, nlen, aseq, 0 );
1055
1056         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
1057         for( i=0; i<njob; i++ ) pair[i][i] = 1;
1058
1059         if( alg == 'H' )
1060         {
1061                 fprintf( stderr, "Calling FOLDALIGN with option '%s'\n", foldalignopt );
1062                 callfoldalign( njob, seq );
1063                 fprintf( stderr, "done.\n" );
1064         }
1065         if( alg == 'B' )
1066         {
1067                 fprintf( stderr, "Running LARA (Bauer et al. http://www.planet-lisa.net/)\n" );
1068                 calllara( njob, seq, "" );
1069                 fprintf( stderr, "done.\n" );
1070         }
1071         if( alg == 'T' )
1072         {
1073                 fprintf( stderr, "Running SLARA (Bauer et al. http://www.planet-lisa.net/)\n" );
1074                 calllara( njob, seq, "-s" );
1075                 fprintf( stderr, "done.\n" );
1076         }
1077         if( alg == 's' )
1078         {
1079                 fprintf( stderr, "Preparing bpp\n" );
1080 //              bpp = AllocateCharCub( njob, nlenmax, 0 );
1081                 bpp = calloc( njob, sizeof( char ** ) );
1082                 preparebpp( njob, bpp );
1083                 fprintf( stderr, "done.\n" );
1084                 fprintf( stderr, "Running MXSCARNA (Tabei et al. http://www.ncrna.org/software/mxscarna)\n" );
1085         }
1086
1087         for( i=0; i<njob; i++ )
1088         {
1089                 pscore = 0.0;
1090                 for( pt=seq[i]; *pt; pt++ )
1091                         pscore += amino_dis[(int)*pt][(int)*pt];
1092                 selfscore[i] = pscore;
1093
1094         }
1095
1096 #if enablemultithread
1097         if( nthread > 0 )
1098         {
1099                 Jobtable jobpos;
1100                 pthread_t *handle;
1101                 pthread_mutex_t mutex_counter;
1102                 pthread_mutex_t mutex_stdout;
1103                 thread_arg_t *targ;
1104
1105                 jobpos.i = 0;
1106                 jobpos.j = 0;
1107
1108                 targ = calloc( nthread, sizeof( thread_arg_t ) );
1109                 handle = calloc( nthread, sizeof( pthread_t ) );
1110                 pthread_mutex_init( &mutex_counter, NULL );
1111                 pthread_mutex_init( &mutex_stdout, NULL );
1112
1113                 for( i=0; i<nthread; i++ )
1114                 {
1115                         targ[i].thread_no = i;
1116                         targ[i].njob = njob;
1117                         targ[i].jobpospt = &jobpos;
1118                         targ[i].name = name;
1119                         targ[i].seq = seq;
1120                         targ[i].localhomtable = localhomtable;
1121                         targ[i].distancemtx = distancemtx;
1122                         targ[i].selfscore = selfscore;
1123                         targ[i].bpp = bpp;
1124                         targ[i].alloclen = alloclen;
1125                         targ[i].mutex_counter = &mutex_counter;
1126                         targ[i].mutex_stdout = &mutex_stdout;
1127
1128 //                      athread( (void *)targ );
1129                         pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
1130 //                      pthread_create( handle+i, NULL, bthread, (void *)(targ+i) );
1131                 }
1132
1133
1134                 for( i=0; i<nthread; i++ )
1135                 {
1136                         pthread_join( handle[i], NULL );
1137                 }
1138                 pthread_mutex_destroy( &mutex_counter );
1139                 pthread_mutex_destroy( &mutex_stdout );
1140                 free( handle );
1141                 free( targ );
1142         }
1143         else
1144 #endif
1145         {
1146                 ilim = njob - 1;
1147                 for( i=0; i<ilim; i++ )
1148                 {
1149                         if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, i+1, 0.0 );
1150                         fprintf( stderr, "% 5d / %d\r", i, njob );
1151                         fflush( stderr );
1152                         for( j=i+1; j<njob; j++ )
1153                         {
1154
1155                                 if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
1156                                 {
1157                                         if( store_dist ) distancemtx[i][j] = 2.0;
1158                                         if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 2.0 );
1159                                         continue;
1160                                 }
1161
1162                                 strcpy( aseq[i], seq[i] );
1163                                 strcpy( aseq[j], seq[j] );
1164                                 clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
1165                                 clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
1166         //                      fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
1167         //                      fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
1168
1169 #if 0
1170                                 fprintf( stderr, "group1 = %.66s", indication1 );
1171                                 fprintf( stderr, "\n" );
1172                                 fprintf( stderr, "group2 = %.66s", indication2 );
1173                                 fprintf( stderr, "\n" );
1174 #endif
1175         //                      for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
1176
1177                                 if( use_fft )
1178                                 {
1179                                         pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
1180 //                                      fprintf( stderr, "pscore (fft) = %f\n", pscore );
1181                                         off1 = off2 = 0;
1182                                 }
1183                                 else
1184                                 {
1185                                         switch( alg )
1186                                         {
1187                                                 case( 'a' ):
1188                                                         pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
1189                                                         off1 = off2 = 0;
1190                                                         break;
1191                                                 case( 't' ):
1192                                                         pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
1193                                                         off1 = off2 = 0;
1194                                                         break;
1195                                                 case( 'A' ):
1196                                                         pscore = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
1197                                                         off1 = off2 = 0;
1198                                                         break;
1199                                                 case( 'N' ):
1200                                                         pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
1201                                                         genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
1202                                                         break;
1203                                                 case( 'K' ):
1204                                                         pscore = genG__align11( mseq1, mseq2, alloclen );
1205                                                         off1 = off2 = 0;
1206                                                         break;
1207                                                 case( 'L' ):
1208                                                         pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
1209                                                         L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
1210                                                         break;
1211                                                 case( 'H' ):
1212                                                         pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
1213                                                         break;
1214                                                 case( 'B' ):
1215                                                 case( 'T' ):
1216                                                         pscore = recalllara( mseq1, mseq2, alloclen );
1217                                                         off1 = off2 = 0;
1218                                                         break;
1219                                                 case( 's' ):
1220                                                         pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
1221                                                         off1 = off2 = 0;
1222                                                         break;
1223                                                 case( 'M' ):
1224                                                         pscore = MSalign11( mseq1, mseq2, alloclen );
1225                                                         break;
1226                                                 default:
1227                                                         ErrorExit( "ERROR IN SOURCE FILE" );
1228                                         }
1229                                 }
1230
1231                                 if( alg == 't' || ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) ) // 't' no jouken ha iranai to omou. if( ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) )
1232                                 {
1233 #if SCOREOUT
1234                                         fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
1235 #endif
1236                                         if( !store_localhom )
1237                                                 ;
1238                                         else if( alg == 'H' )
1239                                                 putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
1240                                         else if( alg != 'S' && alg != 'V' )
1241                                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
1242
1243
1244                                         if( (bunbo=MIN( selfscore[i], selfscore[j] )) == 0.0 || bunbo < pscore )
1245                                                 pscore = 2.0;
1246                                         else
1247                                                 pscore = ( 1.0 - pscore / bunbo ) * 2.0;
1248                                 }
1249                                 else
1250                                 {
1251                                         pscore = 2.0;
1252                                 }
1253
1254                                 if( stdout_align )
1255                                 {
1256                                         if( alg != 't' )
1257                                         {
1258                                                 fprintf( stdout, "sequence %d - sequence %d, pairwise distance = %10.5f\n", i+1, j+1, pscore );
1259                                                 fprintf( stdout, ">%s\n", name[i] );
1260                                                 write1seq( stdout, mseq1[0] );
1261                                                 fprintf( stdout, ">%s\n", name[j] );
1262                                                 write1seq( stdout, mseq2[0] );
1263                                                 fprintf( stdout, "\n" );
1264                                         }
1265                                 }
1266                                 if( stdout_dist ) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, pscore );
1267                                 if( store_dist) distancemtx[i][j] = pscore;
1268                         }
1269                 }
1270         }
1271
1272
1273         if( store_dist )
1274         {
1275                 hat2p = fopen( hat2file, "w" );
1276                 if( !hat2p ) ErrorExit( "Cannot open hat2." );
1277                 WriteHat2_pointer( hat2p, njob, name, distancemtx );
1278                 fclose( hat2p );
1279         }
1280
1281         hat3p = fopen( "hat3", "w" );
1282         if( !hat3p ) ErrorExit( "Cannot open hat3." );
1283         if( store_localhom )
1284         {
1285                 fprintf( stderr, "\n\n##### writing hat3\n" );
1286                 ilim = njob-1;
1287                 for( i=0; i<ilim; i++ )
1288                 {
1289                         for( j=i+1; j<njob; j++ )
1290                         {
1291                                 for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
1292                                 {
1293                                         if( tmpptr->opt == -1.0 ) continue;
1294 // tmptmptmptmptmp
1295 //                                      if( alg == 'B' || alg == 'T' )
1296 //                                              fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, 1.0, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next );
1297 //                                      else
1298                                                 fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2 );
1299                                 }
1300                         }
1301                 }
1302 #if DEBUG
1303                 fprintf( stderr, "calling FreeLocalHomTable\n" );
1304 #endif
1305                 FreeLocalHomTable( localhomtable, njob );
1306 #if DEBUG
1307                 fprintf( stderr, "done. FreeLocalHomTable\n" );
1308 #endif
1309         }
1310         fclose( hat3p );
1311
1312         if( alg == 's' )
1313         {
1314                 char **ptpt;
1315                 for( i=0; i<njob; i++ )
1316                 {
1317                         ptpt = bpp[i];
1318                         while( 1 )
1319                         {
1320                                 if( *ptpt ) free( *ptpt );
1321                                 else break;
1322                                 ptpt++;
1323                         }
1324                         free( bpp[i] );
1325                 }
1326                 free( bpp );
1327         }
1328         free( selfscore );
1329         free( effarr1 );
1330         free( effarr2 );
1331         free( indication1 );
1332         free( indication2 );
1333         if( store_dist ) FreeDoubleMtx( distancemtx );
1334 }
1335
1336 static void WriteOptions( FILE *fp )
1337 {
1338
1339         if( dorp == 'd' ) fprintf( fp, "DNA\n" );
1340         else
1341         {
1342                 if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
1343                 else if( scoremtx ==  1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
1344                 else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
1345         }
1346     fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1347     if( use_fft ) fprintf( fp, "FFT on\n" );
1348
1349         fprintf( fp, "tree-base method\n" );
1350         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
1351         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
1352         if( tbitr || tbweight )
1353         {
1354                 fprintf( fp, "iterate at each step\n" );
1355                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
1356                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
1357                 if( tbweight ) fprintf( fp, "  weighted\n" );
1358                 fprintf( fp, "\n" );
1359         }
1360
1361          fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1362
1363         if( alg == 'a' )
1364                 fprintf( fp, "Algorithm A\n" );
1365         else if( alg == 'A' )
1366                 fprintf( fp, "Algorithm A+\n" );
1367         else if( alg == 'S' )
1368                 fprintf( fp, "Apgorithm S\n" );
1369         else
1370                 fprintf( fp, "Unknown algorithm\n" );
1371
1372     if( use_fft )
1373     {
1374         fprintf( fp, "FFT on\n" );
1375         if( dorp == 'd' )
1376             fprintf( fp, "Basis : 4 nucleotides\n" );
1377         else
1378         {
1379             if( fftscore )
1380                 fprintf( fp, "Basis : Polarity and Volume\n" );
1381             else
1382                 fprintf( fp, "Basis : 20 amino acids\n" );
1383         }
1384         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
1385         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
1386     }
1387         else
1388         fprintf( fp, "FFT off\n" );
1389         fflush( fp );
1390 }
1391
1392
1393 int main( int argc, char *argv[] )
1394 {
1395         int  nlen[M];
1396         char **name, **seq;
1397         char **mseq1, **mseq2;
1398         char **aseq;
1399         char **bseq;
1400         double *eff;
1401         int i;
1402         FILE *infp;
1403         char c;
1404         int alloclen;
1405
1406         arguments( argc, argv );
1407 #ifndef enablemultithread
1408         nthread = 0;
1409 #endif
1410
1411         if( inputfile )
1412         {
1413                 infp = fopen( inputfile, "r" );
1414                 if( !infp )
1415                 {
1416                         fprintf( stderr, "Cannot open %s\n", inputfile );
1417                         exit( 1 );
1418                 }
1419         }
1420         else
1421                 infp = stdin;
1422
1423         getnumlen( infp );
1424         rewind( infp );
1425
1426         if( njob < 2 )
1427         {
1428                 fprintf( stderr, "At least 2 sequences should be input!\n"
1429                                                  "Only %d sequence found.\n", njob );
1430                 exit( 1 );
1431         }
1432         if( njob > M )
1433         {
1434                 fprintf( stderr, "The number of sequences must be < %d\n", M );
1435                 fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
1436                 exit( 1 );
1437         }
1438
1439         alloclen = nlenmax*2;
1440         seq = AllocateCharMtx( njob, alloclen+10 );
1441         aseq = AllocateCharMtx( njob, alloclen+10 );
1442         bseq = AllocateCharMtx( njob, alloclen+10 );
1443         mseq1 = AllocateCharMtx( njob, 0 );
1444         mseq2 = AllocateCharMtx( njob, 0 );
1445         name = AllocateCharMtx( njob, B );
1446
1447         eff = AllocateDoubleVec( njob );
1448
1449 #if 0
1450         Read( name, nlen, seq );
1451 #else
1452         readData_pointer( infp, name, nlen, seq );
1453 #endif
1454         fclose( infp );
1455
1456         constants( njob, seq );
1457
1458 #if 0
1459         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
1460 #endif
1461
1462         initSignalSM();
1463
1464         initFiles();
1465
1466         WriteOptions( trap_g );
1467
1468         c = seqcheck( seq );
1469         if( c )
1470         {
1471                 fprintf( stderr, "Illegal character %c\n", c );
1472                 exit( 1 );
1473         }
1474
1475 //      writePre( njob, name, nlen, seq, 0 );
1476
1477         for( i=0; i<njob; i++ ) eff[i] = 1.0;
1478
1479
1480         for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1481
1482         pairalign( name, nlen, bseq, aseq, mseq1, mseq2, eff, alloclen );
1483
1484         fprintf( trap_g, "done.\n" );
1485 #if DEBUG
1486         fprintf( stderr, "closing trap_g\n" );
1487 #endif
1488         fclose( trap_g );
1489
1490 //      writePre( njob, name, nlen, aseq, !contin );
1491 #if 0
1492         writeData( stdout, njob, name, nlen, aseq );
1493 #endif
1494 #if IODEBUG
1495         fprintf( stderr, "OSHIMAI\n" );
1496 #endif
1497         SHOWVERSION;
1498
1499         if( stdout_dist && nthread > 1 )
1500         {
1501                 fprintf( stderr, "\nThe order of distances is not identical to that in the input file, because of the parallel calculation.  Reorder them by yourself, using sort -n -k 2 | sort -n -k 1 -s\n" );
1502         }
1503         if( stdout_align && nthread > 1 )
1504         {
1505                 fprintf( stderr, "\nThe order of pairwise alignments is not identical to that in the input file, because of the parallel calculation.  Reorder them by yourself.\n" );
1506         }
1507         FreeCharMtx( seq );
1508         FreeCharMtx( aseq );
1509         FreeCharMtx( bseq );
1510         FreeCharMtx( name );
1511         free( mseq1 );
1512         free( mseq2 );
1513         free( eff );
1514
1515         return( 0 );
1516 }