binaries/src/mafft/core/pairlara.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 0
   6
   7 static char *whereispairalign;
   8 static char *laraparams;
   9 static char foldalignopt[1000];
  10
  11 static void t2u( char *seq )
  12 {
  13         while( *seq )
  14         {
  15                 if     ( *seq == 'A' ) *seq = 'a';
  16                 else if( *seq == 'a' ) *seq = 'a';
  17                 else if( *seq == 'T' ) *seq = 'u';
  18                 else if( *seq == 't' ) *seq = 'u';
  19                 else if( *seq == 'U' ) *seq = 'u';
  20                 else if( *seq == 'u' ) *seq = 'u';
  21                 else if( *seq == 'G' ) *seq = 'g';
  22                 else if( *seq == 'g' ) *seq = 'g';
  23                 else if( *seq == 'C' ) *seq = 'c';
  24                 else if( *seq == 'c' ) *seq = 'c';
  25                 else *seq = 'n';
  26                 seq++;
  27         }
  28 }
  29
  30 static float recallpairfoldalign( char **mseq1, char **mseq2, int m1, int m2, int *of1pt, int *of2pt, int alloclen )
  31 {
  32         static FILE *fp = NULL;
  33         float value;
  34         char *aln1;
  35         char *aln2;
  36         int of1tmp, of2tmp;
  37
  38         if( fp == NULL )
  39         {
  40                 fp = fopen( "_foldalignout", "r" );
  41                 if( fp == NULL )
  42                 {
  43                         fprintf( stderr, "Cannot open _foldalignout\n" );
  44                         exit( 1 );
  45                 }
  46         }
  47
  48         aln1 = calloc( alloclen, sizeof( char ) );
  49         aln2 = calloc( alloclen, sizeof( char ) );
  50
  51         readpairfoldalign( fp, *mseq1, *mseq2, aln1, aln2, m1, m2, &of1tmp, &of2tmp, alloclen );
  52
  53         if( strstr( foldalignopt, "-global") )
  54         {
  55                 fprintf( stderr, "Calling G__align11\n" );
  56                 value = G__align11( mseq1, mseq2, alloclen );
  57                 *of1pt = 0;
  58                 *of2pt = 0;
  59         }
  60         else
  61         {
  62                 fprintf( stderr, "Calling L__align11\n" );
  63                 value = L__align11( mseq1, mseq2, alloclen, of1pt, of2pt );
  64         }
  65
  66 //      value = (float)naivepairscore11( *mseq1, *mseq2, penalty ); // nennnotame
  67
  68         if( aln1[0] == 0 )
  69         {
  70                 fprintf( stderr, "FOLDALIGN returned no alignment between %d and %d.  Sequence alignment is used instead.\n", m1+1, m2+1 );
  71         }
  72         else
  73         {
  74                 strcpy( *mseq1, aln1 );
  75                 strcpy( *mseq2, aln2 );
  76                 *of1pt = of1tmp;
  77                 *of2pt = of2tmp;
  78         }
  79
  80 //      value = naivepairscore11( *mseq1, *mseq2, penalty ); // v6.511 ha kore wo tsukau, global nomi dakara.
  81
  82 //      fclose( fp ); // saigo dake yatta houga yoi.
  83
  84 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
  85 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
  86
  87
  88         free( aln1 );
  89         free( aln2 );
  90
  91         return( value );
  92 }
  93
  94 static void callfoldalign( int nseq, char **mseq )
  95 {
  96         FILE *fp;
  97         int i;
  98         int res;
  99         static char com[10000];
 100
 101         for( i=0; i<nseq; i++ )
 102                 t2u( mseq[i] );
 103
 104         fp = fopen( "_foldalignin", "w" );
 105         if( !fp )
 106         {
 107                 fprintf( stderr, "Cannot open _foldalignin\n" );
 108                 exit( 1 );
 109         }
 110         for( i=0; i<nseq; i++ )
 111         {
 112                 fprintf( fp, ">%d\n", i+1 );
 113                 fprintf( fp, "%s\n", mseq[i] );
 114         }
 115         fclose( fp );
 116
 117         sprintf( com, "env PATH=%s  foldalign210 %s _foldalignin > _foldalignout ", whereispairalign, foldalignopt );
 118         res = system( com );
 119         if( res )
 120         {
 121                 fprintf( stderr, "Error in foldalign\n" );
 122                 exit( 1 );
 123         }
 124
 125 }
 126
 127 static void calllara( int nseq, char **mseq, char *laraarg )
 128 {
 129         FILE *fp;
 130         int i;
 131         int res;
 132         static char com[10000];
 133
 134         for( i=0; i<nseq; i++ )
 135
 136         fp = fopen( "_larain", "w" );
 137         if( !fp )
 138         {
 139                 fprintf( stderr, "Cannot open _larain\n" );
 140                 exit( 1 );
 141         }
 142         for( i=0; i<nseq; i++ )
 143         {
 144                 fprintf( fp, ">%d\n", i+1 );
 145                 fprintf( fp, "%s\n", mseq[i] );
 146         }
 147         fclose( fp );
 148
 149
 150 //      fprintf( stderr, "calling LaRA\n" );
 151         sprintf( com, "env PATH=%s:/bin:/usr/bin mafft_lara -i _larain -w _laraout -o _lara.params %s", whereispairalign, laraarg );
 152         res = system( com );
 153         if( res )
 154         {
 155                 fprintf( stderr, "Error in lara\n" );
 156                 exit( 1 );
 157         }
 158 }
 159
 160 static float recalllara( char **mseq1, char **mseq2, int alloclen )
 161 {
 162         static FILE *fp = NULL;
 163         static char *ungap1;
 164         static char *ungap2;
 165         static char *ori1;
 166         static char *ori2;
 167         int res;
 168         static char com[10000];
 169         float value;
 170
 171
 172         if( fp == NULL )
 173         {
 174                 fp = fopen( "_laraout", "r" );
 175                 if( fp == NULL )
 176                 {
 177                         fprintf( stderr, "Cannot open _laraout\n" );
 178                         exit( 1 );
 179                 }
 180                 ungap1 = AllocateCharVec( alloclen );
 181                 ungap2 = AllocateCharVec( alloclen );
 182                 ori1 = AllocateCharVec( alloclen );
 183                 ori2 = AllocateCharVec( alloclen );
 184         }
 185
 186
 187         strcpy( ori1, *mseq1 );
 188         strcpy( ori2, *mseq2 );
 189
 190         fgets( com, 999, fp );
 191         myfgets( com, 9999, fp );
 192         strcpy( *mseq1, com );
 193         myfgets( com, 9999, fp );
 194         strcpy( *mseq2, com );
 195
 196         gappick0( ungap1, *mseq1 );
 197         gappick0( ungap2, *mseq2 );
 198         t2u( ungap1 );
 199         t2u( ungap2 );
 200
 201         if( strcmp( ungap1, ori1 ) || strcmp( ungap2, ori2 ) )
 202         {
 203                 fprintf( stderr, "SEQUENCE CHANGED!!\n" );
 204                 fprintf( stderr, "*mseq1  = %s\n", *mseq1 );
 205                 fprintf( stderr, "ungap1  = %s\n", ungap1 );
 206                 fprintf( stderr, "ori1    = %s\n", ori1 );
 207                 fprintf( stderr, "*mseq2  = %s\n", *mseq2 );
 208                 fprintf( stderr, "ungap2  = %s\n", ungap2 );
 209                 fprintf( stderr, "ori2    = %s\n", ori2 );
 210                 exit( 1 );
 211         }
 212
 213         value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
 214
 215 //      fclose( fp ); // saigo dake yatta houga yoi.
 216
 217         return( value );
 218 }
 219
 220 static float callmxscarna( char **mseq1, char **mseq2, int alloclen )
 221 {
 222         FILE *fp;
 223         int res;
 224         static char com[10000];
 225         float value;
 226
 227
 228         t2u( *mseq1 );
 229         t2u( *mseq2 );
 230         fp = fopen( "_mxscarnain", "w" );
 231         if( !fp )
 232         {
 233                 fprintf( stderr, "Cannot open _mxscarnain\n" );
 234                 exit( 1 );
 235         }
 236         fprintf( fp, ">1\n" );
 237         fprintf( fp, "%s\n", *mseq1 );
 238         fprintf( fp, ">2\n" );
 239         fprintf( fp, "%s\n", *mseq2 );
 240         fclose( fp );
 241
 242         sprintf( com, "env PATH=%s mxscarna _mxscarnain > _mxscarnaout 2>/dev/null", whereispairalign );
 243         res = system( com );
 244         if( res )
 245         {
 246                 fprintf( stderr, "Error in mxscarna\n" );
 247                 exit( 1 );
 248         }
 249
 250         fp = fopen( "_mxscarnaout", "r" );
 251         if( !fp )
 252         {
 253                 fprintf( stderr, "Cannot open _mxscarnaout\n" );
 254                 exit( 1 );
 255         }
 256
 257         fgets( com, 999, fp );
 258         load1SeqWithoutName_new( fp, *mseq1 );
 259         fgets( com, 999, fp );
 260         load1SeqWithoutName_new( fp, *mseq2 );
 261
 262         fclose( fp );
 263
 264 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 265 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 266
 267         value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
 268
 269         return( value );
 270 }
 271
 272
 273 void arguments( int argc, char *argv[] )
 274 {
 275     int c;
 276
 277         foldalignopt[0] = 0;
 278         laraparams = NULL;
 279         inputfile = NULL;
 280         fftkeika = 0;
 281         pslocal = -1000.0;
 282         constraint = 0;
 283         nblosum = 62;
 284         fmodel = 0;
 285         calledByXced = 0;
 286         devide = 0;
 287         use_fft = 0;
 288         fftscore = 1;
 289         fftRepeatStop = 0;
 290         fftNoAnchStop = 0;
 291     weight = 3;
 292     utree = 1;
 293         tbutree = 1;
 294     refine = 0;
 295     check = 1;
 296     cut = 0.0;
 297     disp = 0;
 298     outgap = 1;
 299     alg = 'A';
 300     mix = 0;
 301         tbitr = 0;
 302         scmtd = 5;
 303         tbweight = 0;
 304         tbrweight = 3;
 305         checkC = 0;
 306         treemethod = 'x';
 307         contin = 0;
 308         scoremtx = 1;
 309         kobetsubunkatsu = 0;
 310         divpairscore = 0;
 311         dorp = NOTSPECIFIED;
 312         ppenalty = NOTSPECIFIED;
 313         ppenalty_OP = NOTSPECIFIED;
 314         ppenalty_ex = NOTSPECIFIED;
 315         ppenalty_EX = NOTSPECIFIED;
 316         poffset = NOTSPECIFIED;
 317         kimuraR = NOTSPECIFIED;
 318         pamN = NOTSPECIFIED;
 319         geta2 = GETA2;
 320         fftWinSize = NOTSPECIFIED;
 321         fftThreshold = NOTSPECIFIED;
 322         RNAppenalty = NOTSPECIFIED;
 323         RNApthr = NOTSPECIFIED;
 324
 325     while( --argc > 0 && (*++argv)[0] == '-' )
 326         {
 327         while ( ( c = *++argv[0] ) )
 328                 {
 329             switch( c )
 330             {
 331                                 case 'i':
 332                                         inputfile = *++argv;
 333                                         fprintf( stderr, "inputfile = %s\n", inputfile );
 334                                         --argc;
 335                                         goto nextoption;
 336                                 case 'f':
 337                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 338                                         --argc;
 339                                         goto nextoption;
 340                                 case 'g':
 341                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
 342                                         --argc;
 343                                         goto nextoption;
 344                                 case 'O':
 345                                         ppenalty_OP = (int)( atof( *++argv ) * 1000 - 0.5 );
 346                                         --argc;
 347                                         goto nextoption;
 348                                 case 'E':
 349                                         ppenalty_EX = (int)( atof( *++argv ) * 1000 - 0.5 );
 350                                         --argc;
 351                                         goto nextoption;
 352                                 case 'h':
 353                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
 354                                         --argc;
 355                                         goto nextoption;
 356                                 case 'k':
 357                                         kimuraR = atoi( *++argv );
 358 //                                      fprintf( stderr, "kimuraR = %d\n", kimuraR );
 359                                         --argc;
 360                                         goto nextoption;
 361                                 case 'b':
 362                                         nblosum = atoi( *++argv );
 363                                         scoremtx = 1;
 364 //                                      fprintf( stderr, "blosum %d\n", nblosum );
 365                                         --argc;
 366                                         goto nextoption;
 367                                 case 'j':
 368                                         pamN = atoi( *++argv );
 369                                         scoremtx = 0;
 370                                         TMorJTT = JTT;
 371                                         fprintf( stderr, "jtt %d\n", pamN );
 372                                         --argc;
 373                                         goto nextoption;
 374                                 case 'm':
 375                                         pamN = atoi( *++argv );
 376                                         scoremtx = 0;
 377                                         TMorJTT = TM;
 378                                         fprintf( stderr, "TM %d\n", pamN );
 379                                         --argc;
 380                                         goto nextoption;
 381                                 case 'l':
 382                                         ppslocal = (int)( atof( *++argv ) * 1000 + 0.5 );
 383                                         pslocal = (int)( 600.0 / 1000.0 * ppslocal + 0.5);
 384 //                                      fprintf( stderr, "ppslocal = %d\n", ppslocal );
 385 //                                      fprintf( stderr, "pslocal = %d\n", pslocal );
 386                                         --argc;
 387                                         goto nextoption;
 388                                 case 'd':
 389                                         whereispairalign = *++argv;
 390                                         fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
 391                                         --argc;
 392                                         goto nextoption;
 393                                 case 'p':
 394                                         laraparams = *++argv;
 395                                         fprintf( stderr, "laraparams = %s\n", laraparams );
 396                                         --argc;
 397                                         goto nextoption;
 398 #if 1
 399                                 case 'a':
 400                                         fmodel = 1;
 401                                         break;
 402 #endif
 403                                 case 'r':
 404                                         fmodel = -1;
 405                                         break;
 406                                 case 'D':
 407                                         dorp = 'd';
 408                                         break;
 409                                 case 'P':
 410                                         dorp = 'p';
 411                                         break;
 412                                 case 'e':
 413                                         fftscore = 0;
 414                                         break;
 415 #if 0
 416                                 case 'O':
 417                                         fftNoAnchStop = 1;
 418                                         break;
 419 #endif
 420                                 case 'Q':
 421                                         calledByXced = 1;
 422                                         break;
 423                                 case 'x':
 424                                         disp = 1;
 425                                         break;
 426 #if 0
 427                                 case 'a':
 428                                         alg = 'a';
 429                                         break;
 430 #endif
 431                                 case 'S':
 432                                         alg = 'S';
 433                                         break;
 434                                 case 'L':
 435                                         alg = 'L';
 436                                         break;
 437                                 case 's':
 438                                         alg = 's';
 439                                         break;
 440                                 case 'B':
 441                                         alg = 'B';
 442                                         break;
 443                                 case 'T':
 444                                         alg = 'T';
 445                                         break;
 446                                 case 'H':
 447                                         alg = 'H';
 448                                         break;
 449                                 case 'M':
 450                                         alg = 'M';
 451                                         break;
 452                                 case 'R':
 453                                         alg = 'R';
 454                                         break;
 455                                 case 'N':
 456                                         alg = 'N';
 457                                         break;
 458                                 case 'A':
 459                                         alg = 'A';
 460                                         break;
 461                                 case 'V':
 462                                         alg = 'V';
 463                                         break;
 464                                 case 'C':
 465                                         alg = 'C';
 466                                         break;
 467                                 case 'F':
 468                                         use_fft = 1;
 469                                         break;
 470                                 case 'v':
 471                                         tbrweight = 3;
 472                                         break;
 473                                 case 'y':
 474                                         divpairscore = 1;
 475                                         break;
 476 /* Modified 01/08/27, default: user tree */
 477                                 case 'J':
 478                                         tbutree = 0;
 479                                         break;
 480 /* modification end. */
 481                                 case 'o':
 482 //                                      foldalignopt = *++argv;
 483                                         strcat( foldalignopt, " " );
 484                                         strcat( foldalignopt, *++argv );
 485                                         fprintf( stderr, "foldalignopt = %s\n", foldalignopt );
 486                                         --argc;
 487                                         goto nextoption;
 488                                 case 'z':
 489                                         fftThreshold = atoi( *++argv );
 490                                         --argc;
 491                                         goto nextoption;
 492                                 case 'w':
 493                                         fftWinSize = atoi( *++argv );
 494                                         --argc;
 495                                         goto nextoption;
 496                                 case 'Z':
 497                                         checkC = 1;
 498                                         break;
 499                 default:
 500                     fprintf( stderr, "illegal option %c\n", c );
 501                     argc = 0;
 502                     break;
 503             }
 504                 }
 505                 nextoption:
 506                         ;
 507         }
 508     if( argc == 1 )
 509     {
 510         cut = atof( (*argv) );
 511         argc--;
 512     }
 513     if( argc != 0 )
 514     {
 515         fprintf( stderr, "options: Check source file !\n" );
 516         exit( 1 );
 517     }
 518         if( tbitr == 1 && outgap == 0 )
 519         {
 520                 fprintf( stderr, "conflicting options : o, m or u\n" );
 521                 exit( 1 );
 522         }
 523         if( alg == 'C' && outgap == 0 )
 524         {
 525                 fprintf( stderr, "conflicting options : C, o\n" );
 526                 exit( 1 );
 527         }
 528 }
 529
 530 int countamino( char *s, int end )
 531 {
 532         int val = 0;
 533         while( end-- )
 534                 if( *s++ != '-' ) val++;
 535         return( val );
 536 }
 537
 538 static void pairalign( char name[M][B], int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, double *effarr, int alloclen )
 539 {
 540         int i, j, ilim;
 541         int clus1, clus2;
 542         int off1, off2;
 543         float pscore = 0.0; // by D.Mathog
 544         static char *indication1, *indication2;
 545         FILE *hat2p, *hat3p;
 546         static double **distancemtx;
 547         static double *effarr1 = NULL;
 548         static double *effarr2 = NULL;
 549         char *pt;
 550         char *hat2file = "hat2";
 551         LocalHom **localhomtable, *tmpptr;
 552         static char **pair;
 553         int intdum;
 554         double bunbo;
 555
 556         localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 557         for( i=0; i<njob; i++)
 558         {
 559                 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
 560                 for( j=0; j<njob; j++)
 561                 {
 562                         localhomtable[i][j].start1 = -1;
 563                         localhomtable[i][j].end1 = -1;
 564                         localhomtable[i][j].start2 = -1;
 565                         localhomtable[i][j].end2 = -1;
 566                         localhomtable[i][j].opt = -1.0;
 567                         localhomtable[i][j].next = NULL;
 568                         localhomtable[i][j].nokori = 0;
 569                 }
 570         }
 571
 572         if( effarr1 == NULL )
 573         {
 574                 distancemtx = AllocateDoubleMtx( njob, njob );
 575                 effarr1 = AllocateDoubleVec( njob );
 576                 effarr2 = AllocateDoubleVec( njob );
 577                 indication1 = AllocateCharVec( 150 );
 578                 indication2 = AllocateCharVec( 150 );
 579 #if 0
 580 #else
 581                 pair = AllocateCharMtx( njob, njob );
 582 #endif
 583         }
 584
 585 #if 0
 586         fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
 587 #endif
 588
 589 #if 0
 590         for( i=0; i<njob; i++ )
 591                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
 592 #endif
 593
 594
 595 //      writePre( njob, name, nlen, aseq, 0 );
 596
 597         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
 598         for( i=0; i<njob; i++ ) pair[i][i] = 1;
 599
 600         if( alg == 'H' )
 601         {
 602                 fprintf( stderr, "Calling FOLDALIGN with option '%s'\n", foldalignopt );
 603                 callfoldalign( njob, seq );
 604                 fprintf( stderr, "done.\n" );
 605         }
 606         if( alg == 'B' )
 607         {
 608                 fprintf( stderr, "Calling LARA\n" );
 609                 calllara( njob, seq, "" );
 610                 fprintf( stderr, "done.\n" );
 611         }
 612         if( alg == 'T' )
 613         {
 614                 fprintf( stderr, "Calling SLARA\n" );
 615                 calllara( njob, seq, "-s" );
 616                 fprintf( stderr, "done.\n" );
 617         }
 618
 619         ilim = njob - 1;
 620         for( i=0; i<ilim; i++ )
 621         {
 622                 fprintf( stderr, "% 5d / %d\r", i, njob );
 623                 for( j=i+1; j<njob; j++ )
 624                 {
 625
 626                         if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
 627                         {
 628                                 distancemtx[i][j] = pscore;
 629                                 continue;
 630                         }
 631
 632                         strcpy( aseq[i], seq[i] );
 633                         strcpy( aseq[j], seq[j] );
 634                         clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
 635                         clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
 636         //              fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
 637         //              fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
 638
 639 #if 0
 640                         fprintf( stderr, "group1 = %.66s", indication1 );
 641                         fprintf( stderr, "\n" );
 642                         fprintf( stderr, "group2 = %.66s", indication2 );
 643                         fprintf( stderr, "\n" );
 644 #endif
 645         //              for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
 646
 647 #if 1
 648                         if( use_fft )
 649                         {
 650                                 pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
 651                                 off1 = off2 = 0;
 652                         }
 653                         else
 654 #endif
 655                         {
 656                                 switch( alg )
 657                                 {
 658                                         case( 'a' ):
 659                                                 pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 660                                                 off1 = off2 = 0;
 661                                                 break;
 662                                         case( 'A' ):
 663                                                 pscore = G__align11( mseq1, mseq2, alloclen, NULL, 0, NULL );
 664                                                 off1 = off2 = 0;
 665                                                 break;
 666 #if 0
 667                                         case( 'V' ):
 668                                                 pscore = VAalign11( mseq1, mseq2, alloclen, &off1, &off2, localhomtable[i]+j );
 669                                                 fprintf( stderr, "i,j = %d,%d, score = %f\n", i,j, pscore );
 670                                                 break;
 671                                         case( 'S' ):
 672                                                 fprintf( stderr, "aligning %d-%d\n", i, j );
 673                                                 pscore = suboptalign11( mseq1, mseq2, alloclen, &off1, &off2, localhomtable[i]+j );
 674                                                 fprintf( stderr, "i,j = %d,%d, score = %f\n", i,j, pscore );
 675                                                 break;
 676 #endif
 677                                         case( 'N' ):
 678                                                 pscore = genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
 679 //                                              fprintf( stderr, "pscore = %f\n", pscore );
 680                                                 break;
 681                                         case( 'L' ):
 682                                                 pscore = L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
 683 //                                              fprintf( stderr, "pscore (1) = %f\n", pscore );
 684 //                                              pscore = (float)naivepairscore11( *mseq1, *mseq2, penalty ); // nennnotame
 685 //                                              fprintf( stderr, "pscore (2) = %f\n\n", pscore );
 686                                                 break;
 687                                         case( 'H' ):
 688                                                 pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
 689                                                 break;
 690                                         case( 'B' ):
 691                                         case( 'T' ):
 692                                                 pscore = recalllara( mseq1, mseq2, alloclen );
 693                                                 off1 = off2 = 0;
 694 //                                              fprintf( stderr, "lara, pscore = %f\n", pscore );
 695                                                 break;
 696                                         case( 's' ):
 697                                                 pscore = callmxscarna( mseq1, mseq2, alloclen );
 698                                                 off1 = off2 = 0;
 699 //                                              fprintf( stderr, "scarna, pscore = %f\n", pscore );
 700                                                 break;
 701                                         case( 'M' ):
 702 //                                              pscore = MSalign11( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL );
 703                                                 pscore = MSalign11( mseq1, mseq2, alloclen );
 704 //                                              fprintf( stderr, "pscore = %f\n", pscore );
 705                                                 break;
 706                                                 ErrorExit( "ERROR IN SOURCE FILE" );
 707                                 }
 708                         }
 709                         distancemtx[i][j] = pscore;
 710 #if SCOREOUT
 711                         fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
 712 #endif
 713 //                      fprintf( stderr, "pslocal = %d\n", pslocal );
 714 //                      offset = makelocal( *mseq1, *mseq2, pslocal );
 715 #if 0
 716                         fprintf( stderr, "off1 = %d, off2 = %d\n", off1, off2 );
 717                         fprintf( stderr, ">%d\n%s\n>%d\n%s\n>\n", i, mseq1[0], j, mseq2[0] );
 718 #endif
 719
 720 //                      putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, countamino( *mseq1, off1 ), countamino( *mseq2, off2 ), pscore, strlen( mseq1[0] ) );
 721 //                      fprintf( stderr, "pscore = %f\n", pscore );
 722                         if( alg == 'H' )
 723 //                      if( alg == 'H' || alg == 's' || alg == 'B' ) // next version
 724                                 putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
 725                         else if( alg != 'S' && alg != 'V' )
 726                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
 727                 }
 728         }
 729         for( i=0; i<njob; i++ )
 730         {
 731                 pscore = 0.0;
 732                 for( pt=seq[i]; *pt; pt++ )
 733                         pscore += amino_dis[(int)*pt][(int)*pt];
 734                 distancemtx[i][i] = pscore;
 735
 736         }
 737
 738         ilim = njob-1;
 739         for( i=0; i<ilim; i++ )
 740         {
 741                 for( j=i+1; j<njob; j++ )
 742                 {
 743                         bunbo = MIN( distancemtx[i][i], distancemtx[j][j] );
 744                         if( bunbo == 0.0 )
 745                                 distancemtx[i][j] = 2.0;
 746                         else
 747                                 distancemtx[i][j] = ( 1.0 - distancemtx[i][j] / bunbo ) * 2.0;
 748                 }
 749         }
 750
 751         hat2p = fopen( hat2file, "w" );
 752         if( !hat2p ) ErrorExit( "Cannot open hat2." );
 753         WriteHat2( hat2p, njob, name, distancemtx );
 754         fclose( hat2p );
 755
 756         fprintf( stderr, "##### writing hat3\n" );
 757         hat3p = fopen( "hat3", "w" );
 758         if( !hat3p ) ErrorExit( "Cannot open hat3." );
 759         ilim = njob-1;
 760         for( i=0; i<ilim; i++ )
 761         {
 762                 for( j=i+1; j<njob; j++ )
 763                 {
 764                         for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
 765                         {
 766                                 if( tmpptr->opt == -1.0 ) continue;
 767                                 fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next );
 768                         }
 769                 }
 770         }
 771         fclose( hat3p );
 772 #if DEBUG
 773         fprintf( stderr, "calling FreeLocalHomTable\n" );
 774 #endif
 775         FreeLocalHomTable( localhomtable, njob );
 776 #if DEBUG
 777         fprintf( stderr, "done. FreeLocalHomTable\n" );
 778 #endif
 779 }
 780
 781 static void WriteOptions( FILE *fp )
 782 {
 783
 784         if( dorp == 'd' ) fprintf( fp, "DNA\n" );
 785         else
 786         {
 787                 if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
 788                 else if( scoremtx ==  1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
 789                 else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
 790         }
 791     fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 792     if( use_fft ) fprintf( fp, "FFT on\n" );
 793
 794         fprintf( fp, "tree-base method\n" );
 795         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
 796         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
 797         if( tbitr || tbweight )
 798         {
 799                 fprintf( fp, "iterate at each step\n" );
 800                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
 801                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
 802                 if( tbweight ) fprintf( fp, "  weighted\n" );
 803                 fprintf( fp, "\n" );
 804         }
 805
 806          fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 807
 808         if( alg == 'a' )
 809                 fprintf( fp, "Algorithm A\n" );
 810         else if( alg == 'A' )
 811                 fprintf( fp, "Algorithm A+\n" );
 812         else if( alg == 'S' )
 813                 fprintf( fp, "Apgorithm S\n" );
 814         else if( alg == 'C' )
 815                 fprintf( fp, "Apgorithm A+/C\n" );
 816         else
 817                 fprintf( fp, "Unknown algorithm\n" );
 818
 819     if( use_fft )
 820     {
 821         fprintf( fp, "FFT on\n" );
 822         if( dorp == 'd' )
 823             fprintf( fp, "Basis : 4 nucleotides\n" );
 824         else
 825         {
 826             if( fftscore )
 827                 fprintf( fp, "Basis : Polarity and Volume\n" );
 828             else
 829                 fprintf( fp, "Basis : 20 amino acids\n" );
 830         }
 831         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
 832         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
 833     }
 834         else
 835         fprintf( fp, "FFT off\n" );
 836         fflush( fp );
 837 }
 838
 839
 840 int main( int argc, char *argv[] )
 841 {
 842         static int  nlen[M];
 843         static char name[M][B], **seq;
 844         static char **mseq1, **mseq2;
 845         static char **aseq;
 846         static char **bseq;
 847         static double *eff;
 848         int i;
 849         FILE *infp;
 850         char c;
 851         int alloclen;
 852
 853         arguments( argc, argv );
 854
 855         if( inputfile )
 856         {
 857                 infp = fopen( inputfile, "r" );
 858                 if( !infp )
 859                 {
 860                         fprintf( stderr, "Cannot open %s\n", inputfile );
 861                         exit( 1 );
 862                 }
 863         }
 864         else
 865                 infp = stdin;
 866
 867         getnumlen( infp );
 868         rewind( infp );
 869
 870         if( njob < 2 )
 871         {
 872                 fprintf( stderr, "At least 2 sequences should be input!\n"
 873                                                  "Only %d sequence found.\n", njob );
 874                 exit( 1 );
 875         }
 876         if( njob > M )
 877         {
 878                 fprintf( stderr, "The number of sequences must be < %d\n", M );
 879                 fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
 880                 exit( 1 );
 881         }
 882
 883         seq = AllocateCharMtx( njob, nlenmax*9+1 );
 884         aseq = AllocateCharMtx( njob, nlenmax*9+1 );
 885         bseq = AllocateCharMtx( njob, nlenmax*9+1 );
 886         mseq1 = AllocateCharMtx( njob, 0 );
 887         mseq2 = AllocateCharMtx( njob, 0 );
 888         alloclen = nlenmax*9;
 889
 890         eff = AllocateDoubleVec( njob );
 891
 892 #if 0
 893         Read( name, nlen, seq );
 894 #else
 895         readData( infp, name, nlen, seq );
 896 #endif
 897         fclose( infp );
 898
 899         constants( njob, seq );
 900
 901 #if 0
 902         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
 903 #endif
 904
 905         initSignalSM();
 906
 907         initFiles();
 908
 909         WriteOptions( trap_g );
 910
 911         c = seqcheck( seq );
 912         if( c )
 913         {
 914                 fprintf( stderr, "Illegal character %c\n", c );
 915                 exit( 1 );
 916         }
 917
 918 //      writePre( njob, name, nlen, seq, 0 );
 919
 920         for( i=0; i<njob; i++ ) eff[i] = 1.0;
 921
 922
 923         for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
 924
 925         pairalign( name, nlen, bseq, aseq, mseq1, mseq2, eff, alloclen );
 926
 927         fprintf( trap_g, "done.\n" );
 928 #if DEBUG
 929         fprintf( stderr, "closing trap_g\n" );
 930 #endif
 931         fclose( trap_g );
 932
 933 //      writePre( njob, name, nlen, aseq, !contin );
 934 #if 0
 935         writeData( stdout, njob, name, nlen, aseq );
 936 #endif
 937 #if IODEBUG
 938         fprintf( stderr, "OSHIMAI\n" );
 939 #endif
 940         SHOWVERSION;
 941         return( 0 );
 942 }