binaries/src/mafft/core/pairash.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 0
   6
   7 static int usecache;
   8 static char *whereispairalign;
   9 static double scale;
  10 static int *alreadyoutput;
  11 static int equivthreshold;
  12 static int equivwinsize;
  13 static int equivshortestlen;
  14
  15 static void cutpath( char *s )
  16 {
  17         char *pos;
  18         pos = s + strlen( s );
  19
  20         while( --pos >= s )
  21         {
  22                 if( *pos == '/' ) break;
  23         }
  24
  25         strcpy( s, pos+1 );
  26 }
  27
  28 static char getchainid( char *s )
  29 {
  30         s += strlen( s ) - 2;
  31         if( isspace( s[0] ) && isalnum( s[1] ) )
  32                 return( s[1] );
  33         else
  34                 return( 'A' );
  35 }
  36
  37 static void extractfirstword( char *s )
  38 {
  39         while( *s )
  40         {
  41                 if( isspace( *s ) ) break;
  42                 s++;
  43         }
  44         *s = 0;
  45 }
  46
  47 static char *strip( char *s )
  48 {
  49         char *v;
  50
  51         while( *s )
  52         {
  53                 if( !isspace( *s ) ) break;
  54                 s++;
  55         }
  56         v = s;
  57
  58         s += strlen( v ) - 1;
  59         while( s>=v )
  60         {
  61                 if( !isspace( *s ) )
  62                 {
  63                         *(s+1) = 0;
  64                         break;
  65                 }
  66                 s--;
  67         }
  68
  69         return( v );
  70 }
  71
  72 #if 0
  73 static void makeequivdouble( double *d, char *c )
  74 {
  75         while( *c )
  76         {
  77                 *d++ = (double)( *c++ - '0' );
  78         }
  79 }
  80
  81 static void maskequiv( double *d, int n )
  82 {
  83         int halfwin;
  84         int ok;
  85         int i, j;
  86
  87         halfwin = (int)( equivwinsize / 2 );
  88
  89         for( i=0; i<n; i++ )
  90         {
  91                 ok = 1;
  92                 for( j = i-halfwin; j<i+halfwin; j++ )
  93                 {
  94                         if( j<0 || n=<j ) continue;
  95                         if( d[j] <= 0.0 )
  96                         {
  97                                 ok = 0;
  98                                 break;
  99                         }
 100                 }
 101                 if( ok == 0 ) d[i] = 0.0;
 102         }
 103 }
 104 #else
 105 static void maskequiv( double *d, int n )
 106 {
 107         int i, len;
 108         int count;
 109         len = 0;
 110         double *dbk, *dori, *dbkori;
 111
 112         dbk = calloc( n, sizeof( double ) );
 113
 114         dbkori = dbk;
 115         dori = d;
 116         count = n;
 117         while( count-- )
 118         {
 119                 *dbk++ = *d++;
 120         }
 121
 122         dbk = dbkori;
 123         d = dori;
 124         len = 0;
 125
 126
 127         for( i=0; i<n; i++ )
 128         {
 129                 if( d[i] > 0.0 )
 130                 {
 131                         len += 1;
 132                         d[i] = 0.0;
 133                 }
 134                 else
 135                 {
 136                         d[i] = 0.0;
 137                         if( len >= equivshortestlen )
 138                         {
 139                                 len++;
 140                                 while( len-- ) d[i-len] = dbk[i-len];
 141                         }
 142                         len = 0;
 143                 }
 144         }
 145
 146         if( len >= equivshortestlen )
 147         {
 148                 len++;
 149                 while( len-- ) d[n-len] = dbk[n-len];
 150         }
 151
 152         free( dbk );
 153 }
 154 #endif
 155
 156 static void makeequivdouble_tmalign( double *d, char *c, int n )
 157 {
 158         double tmpd;
 159         double *dbk;
 160         int tmpi;
 161         char s;
 162         dbk = d;
 163         while( *c )
 164         {
 165                 if( ( s=*c++ ) == ':' )
 166                         tmpi = 9;
 167                 else if( s == '.' )
 168                         tmpi = 4;
 169                 else
 170                         tmpi = 0;
 171 //              tmpd = (double)( tmpi + 1 - equivthreshold ) / ( 10 - equivthreshold ) * 9.0;
 172 //              if( tmpd < 0.0 ) tmpd = 0.0;
 173                 tmpd = (double)( tmpi );
 174 //              *d++ = (int)tmpd;
 175                 *d++ = tmpd;
 176         }
 177
 178         d = dbk;
 179 //      maskequiv( d, n );
 180 }
 181
 182 static void makeequivdouble_threshold( double *d, char *c, int n )
 183 {
 184         double tmpd;
 185         double *dbk;
 186         int tmpi;
 187         dbk = d;
 188         while( *c )
 189         {
 190                 tmpi = (int)( *c++ - '0' );
 191                 tmpd = (double)( tmpi + 1 - equivthreshold ) / ( 10 - equivthreshold ) * 9.0;
 192                 if( tmpd < 0.0 ) tmpd = 0.0;
 193 //              *d++ = (int)tmpd;
 194                 *d++ = tmpd;
 195         }
 196
 197         d = dbk;
 198         maskequiv( d, n );
 199 }
 200
 201 static void readtmalign( FILE *fp, char *seq1, char *seq2, double *equiv )
 202 {
 203         static char *line = NULL;
 204         static char *equivchar = NULL;
 205         int n;
 206
 207
 208         if( equivchar == NULL )
 209         {
 210                 equivchar = calloc( nlenmax * 2 + 1, sizeof( char ) );
 211                 line = calloc( nlenmax * 2 + 1, sizeof( char ) );
 212         }
 213         seq1[0] = 0;
 214         seq2[0] = 0;
 215         equivchar[0] = 0;
 216
 217
 218 //      system( "vi _tmalignout" );
 219         while( 1 )
 220         {
 221                 if( feof( fp ) )
 222                 {
 223                         fprintf( stderr, "Error in TMalign\n" );
 224                         exit( 1 );
 225                 }
 226                 fgets( line, 999, fp );
 227 //              fprintf( stdout, "line = :%s:\n", line );
 228                 if( !strncmp( line+5, "denotes the residue pairs", 20 ) ) break;
 229         }
 230         fgets( line, nlenmax*2, fp );
 231         strcat( seq1, strip( line ) );
 232
 233         fgets( line, nlenmax*2, fp );
 234         strcat( equivchar, strip( line ) );
 235
 236         fgets( line, nlenmax*2, fp );
 237         strcat( seq2, strip( line ) );
 238
 239 #if 0
 240         printf( "seq1=%s\n", seq1 );
 241         printf( "seq2=%s\n", seq2 );
 242         printf( "equi=%s\n", equivchar );
 243 exit( 1 );
 244 #endif
 245         n = strlen( seq1 );
 246         makeequivdouble_tmalign( equiv, equivchar, n );
 247
 248 #if 0
 249         fprintf( stdout, "\n" );
 250         for( i=0; i<n; i++ )
 251         {
 252                 fprintf( stdout, "%1.0f", equiv[i] );
 253         }
 254         fprintf( stdout, "\n" );
 255         exit( 1 );
 256 #endif
 257 }
 258 static void readrash( FILE *fp, char *seq1, char *seq2, double *equiv )
 259 {
 260         char line[1000];
 261         static char *equivchar = NULL;
 262         int n;
 263
 264
 265         if( equivchar == NULL )
 266         {
 267                 equivchar = calloc( nlenmax * 2, sizeof( char ) );
 268         }
 269         seq1[0] = 0;
 270         seq2[0] = 0;
 271         equivchar[0] = 0;
 272
 273         while( 1 )
 274         {
 275                 fgets( line, 999, fp );
 276 //              fprintf( stdout, "line = :%s:\n", line );
 277                 if( !strncmp( line, "Query ", 6 ) ) break;
 278                 if( feof( fp ) ) break;
 279                 if( !strncmp( line, "QUERY ", 6 ) )
 280                 {
 281                         strcat( seq1, strip( line+5 ) );
 282                         fgets( line, 999, fp );
 283                 }
 284                 if( !strncmp( line, "TEMPL ", 6 ) )
 285                 {
 286                         strcat( seq2, strip( line+5 ) );
 287                         fgets( line, 999, fp );
 288                 }
 289                 if( !strncmp( line, "Equiva", 6 ) )
 290                 {
 291                         strcat( equivchar, strip( line+11 ) );
 292                         fgets( line, 999, fp );
 293                 }
 294         }
 295 #if 0
 296         printf( "seq1=:%s:\n", seq1 );
 297         printf( "seq2=:%s:\n", seq2 );
 298         printf( "equi=:%s:\n", equivchar );
 299 exit( 1 );
 300 #endif
 301         n = strlen( seq1 );
 302         makeequivdouble_threshold( equiv, equivchar, n );
 303
 304 #if 0
 305         fprintf( stdout, "\n" );
 306         for( i=0; i<n; i++ )
 307         {
 308                 fprintf( stdout, "%1.0f", equiv[i] );
 309         }
 310         fprintf( stdout, "\n" );
 311 #endif
 312 }
 313
 314 static int checkcbeta( FILE *fp )
 315 {
 316         char linec[1000];
 317         while( 1 )
 318         {
 319                 fgets( linec, 999, fp );
 320                 if( feof( fp ) ) break;
 321                 if( !strncmp( "ATOM ", linec, 5 ) )
 322                 {
 323                         if( !strncmp( "CB ", linec+13, 3 ) ) return( 0 );
 324                 }
 325         }
 326         return( 1 );
 327 }
 328
 329
 330 static float calltmalign( char **mseq1, char **mseq2, double *equiv, char *fname1, char *chain1, char *fname2, char *chain2, int alloclen )
 331 {
 332         FILE *fp;
 333         int res;
 334         static char com[10000];
 335         float value;
 336         char cachedir[10000];
 337         char cachefile[10000];
 338         int runnow;
 339
 340
 341         if( usecache )
 342         {
 343                 sprintf( cachedir, "%s/.tmalignoutcache", getenv( "HOME" ) );
 344                 sprintf( com, "mkdir -p %s", cachedir );
 345                 system( com );
 346
 347                 sprintf( cachefile, "%s/%s%s-%s%s", cachedir, fname1, chain1, fname2, chain2 );
 348
 349                 runnow = 0;
 350                 fp = fopen( cachefile, "r" );
 351                 if( fp == NULL ) runnow = 1;
 352                 else
 353                 {
 354                         fgets( com, 100, fp );
 355                         if( strncmp( com, "successful", 10 ) ) runnow = 1;
 356                         fclose( fp );
 357                 }
 358         }
 359         else
 360         {
 361                 runnow = 1;
 362         }
 363
 364         if( runnow )
 365         {
 366 #if 0
 367                 sprintf( com, "ln -s %s %s.pdb 2>_dum", fname1, fname1 );
 368                 res = system( com );
 369                 sprintf( com, "ln -s %s %s.pdb 2>_dum", fname2, fname2 );
 370                 res = system( com );
 371 #endif
 372                 sprintf( com, "\"%s/TMalign\"  %s.pdb %s.pdb > _tmalignout 2>_dum", whereispairalign, fname1, fname2 );
 373                 fprintf( stderr, "command = %s\n", com );
 374                 res = system( com );
 375                 if( res )
 376                 {
 377                         fprintf( stderr, "Error in TMalign\n" );
 378                         exit( 1 );
 379                 }
 380
 381         }
 382         else
 383         {
 384                 fprintf( stderr, "Cache is not supported!\n" );
 385                 exit( 1 );
 386         }
 387
 388         fp = fopen( "_tmalignout", "r" );
 389         if( !fp )
 390         {
 391                 fprintf( stderr, "Cannot open _tmalignout\n" );
 392                 exit( 1 );
 393         }
 394
 395         readtmalign( fp, *mseq1, *mseq2, equiv );
 396
 397         fclose( fp );
 398
 399 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 400 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 401
 402         value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
 403
 404         return( value );
 405 }
 406
 407 static float callrash( int mem1, int mem2, char **mseq1, char **mseq2, double *equiv, char *fname1, char *chain1, char *fname2, char *chain2, int alloclen )
 408 {
 409         FILE *fp;
 410         int res;
 411         static char com[10000];
 412         float value;
 413         char cachedir[10000];
 414         char cachefile[10000];
 415         int runnow;
 416         char pairid[1000];
 417
 418         sprintf( pairid, "%d-%d", mem1, mem2 );
 419 //      fprintf( stderr, "pairid = %s\n", pairid );
 420
 421         if( usecache )
 422         {
 423                 sprintf( cachedir, "%s/.rashoutcache", getenv( "HOME" ) );
 424                 sprintf( com, "mkdir -p %s", cachedir );
 425                 system( com );
 426
 427                 sprintf( cachefile, "%s/%s%s-%s%s", cachedir, fname1, chain1, fname2, chain2 );
 428
 429                 runnow = 0;
 430                 fp = fopen( cachefile, "r" );
 431                 if( fp == NULL ) runnow = 1;
 432                 else
 433                 {
 434                         fgets( com, 100, fp );
 435                         if( strncmp( com, "successful", 10 ) ) runnow = 1;
 436                         fclose( fp );
 437                 }
 438         }
 439         else
 440         {
 441                 runnow = 1;
 442         }
 443
 444         if( runnow )
 445         {
 446 #if 0
 447                 sprintf( com, "ln -s %s %s.pdb 2>_dum", fname1, fname1 );
 448                 res = system( com );
 449                 sprintf( com, "ln -s %s %s.pdb 2>_dum", fname2, fname2 );
 450                 res = system( com );
 451 #endif
 452 #if 0  // 091127, pdp nai!
 453                 sprintf( com, "env PATH=%s PDP_ASH.pl --qf %s.pdb --qc %s --tf %s.pdb --tc %s > _rashout 2>_dum", whereispairalign, fname1, chain1, fname2, chain2 );
 454 #else
 455                 sprintf( com, "\"%s/rash\" --qf %s.pdb --qc %s --tf %s.pdb --tc %s --of %s.pdbpair > %s.rashout 2>%s.dum", whereispairalign, fname1, chain1, fname2, chain2, pairid, pairid, pairid );
 456 #endif
 457                 fprintf( stderr, "command = %s\n", com );
 458                 res = system( com );
 459                 if( res )
 460                 {
 461                         fprintf( stderr, "Error in structural alignment\n" );
 462                         exit( 1 );
 463                 }
 464                 sprintf( com, "awk '/^REMARK/,/^TER/' %s.pdbpair > %s.%s-x-%s.%s.pdbpair", pairid, fname1, chain1, fname2, chain2 );
 465                 res = system( com );
 466
 467                 sprintf( com, "awk '/^REMARK/,/^TER/{next} 1' %s.pdbpair > %s.%s-x-%s.%s.pdbpair", pairid, fname2, chain2, fname1, chain1 );
 468                 res = system( com );
 469
 470                 sprintf( com, "rm %s.pdbpair", pairid );
 471                 res = system( com );
 472
 473
 474         }
 475         else
 476         {
 477                 fprintf( stderr, "Use cache!\n" );
 478                 sprintf( com, "grep -v successful %s > %s.rashout", cachefile, pairid );
 479                 system( com );
 480         }
 481
 482         if( usecache && runnow )
 483         {
 484                 sprintf( com, "echo successful >  %s", cachefile );
 485                 system( com );
 486                 sprintf( com, "cat %s.rashout >>  %s", pairid, cachefile );
 487                 system( com );
 488         }
 489
 490         sprintf( com, "%s.rashout", pairid );
 491         fp = fopen( com, "r" );
 492         if( !fp )
 493         {
 494                 fprintf( stderr, "Cannot open %s\n", com );
 495                 exit( 1 );
 496         }
 497
 498         readrash( fp, *mseq1, *mseq2, equiv );
 499
 500         fclose( fp );
 501
 502 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 503 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 504
 505
 506         value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
 507
 508         return( value );
 509 }
 510
 511 static void preparetmalign( FILE *fp, char ***strfiles, char ***chainids, char ***seqpt, char ***mseq1pt, char ***mseq2pt, double **equivpt, int *alloclenpt )
 512 {
 513         int i, res;
 514         char *dumseq;
 515         char line[1000];
 516         char fname[1000];
 517         char command[1000];
 518         int linenum, istr, nstr;
 519         FILE *checkfp;
 520         char *sline;
 521         int use[1000];
 522         linenum = 0;
 523         nstr = 0;
 524         while( 1 )
 525         {
 526                 fgets( line, 999, fp );
 527                 if( feof( fp ) ) break;
 528                 sline = strip( line );
 529                 use[linenum] = 1;
 530                 if( sline[0] == '#' || strlen( sline ) < 2 )
 531                 {
 532                         use[linenum] = 0;
 533                         linenum++;
 534                         continue;
 535                 }
 536                 extractfirstword( sline );
 537                 checkfp = fopen( sline, "r" );
 538                 if( checkfp == NULL )
 539                 {
 540                         fprintf( stderr, "Cannot open %s.\n", sline );
 541                         exit( 1 );
 542                 }
 543 #if 0
 544                 fgets( linec, 999, checkfp );
 545                 if( strncmp( "HEADER ", linec, 7 ) )
 546                 {
 547                         fprintf( stderr, "Check the format of %s.\n", sline );
 548                         exit( 1 );
 549                 }
 550 #endif
 551                 if( checkcbeta( checkfp ) )
 552                 {
 553                         fprintf( stderr, "%s has no C-beta atoms.\n", sline );
 554                         exit( 1 );
 555                 }
 556                 else
 557                         nstr++;
 558                 fclose( checkfp );
 559                 linenum++;
 560         }
 561         njob = nstr;
 562         fprintf( stderr, "nstr = %d\n", nstr );
 563
 564         *strfiles = AllocateCharMtx( nstr, 1000 );
 565         *chainids = AllocateCharMtx( nstr, 2 );
 566
 567         rewind( fp );
 568         istr = 0;
 569         linenum = 0;
 570         while( 1 )
 571         {
 572                 fgets( line, 999, fp );
 573                 if( feof( fp ) ) break;
 574                 sline = strip( line );
 575                 if( use[linenum++] )
 576                 {
 577                         (*chainids)[istr][0] = getchainid( sline );
 578                         (*chainids)[istr][1] = 0;
 579                         extractfirstword( sline );
 580                         sprintf( fname, "%s", sline );
 581                         cutpath( fname );
 582                         sprintf( command, "cp %s %s.pdb", sline, fname );
 583                         system( command );
 584                         sprintf( command, "perl \"%s/clean.pl\" %s.pdb", whereispairalign, fname );
 585                         res = system( command );
 586                         if( res )
 587                         {
 588                                 fprintf( stderr, "error: Install clean.pl\n" );
 589                                 exit( 1 );
 590                         }
 591                         strcpy( (*strfiles)[istr++], fname );
 592                 }
 593         }
 594
 595         *seqpt = AllocateCharMtx( njob, nlenmax*2+1 );
 596         *mseq1pt = AllocateCharMtx( njob, 0 );
 597         *mseq2pt = AllocateCharMtx( njob, 0 );
 598         *equivpt = AllocateDoubleVec( nlenmax*2+1 );
 599         *alloclenpt = nlenmax*2;
 600         dumseq = AllocateCharVec( nlenmax*2+1 );
 601         alreadyoutput = AllocateIntVec( njob );
 602         for( i=0; i<njob; i++ ) alreadyoutput[i] = 0;
 603
 604         for( i=0; i<istr; i++ )
 605         {
 606                 fprintf( stderr, "i=%d\n", i );
 607                 (*seqpt)[i][0] = 0;
 608
 609                 (*mseq1pt)[0] = (*seqpt)[i];
 610                 (*mseq2pt)[0] = dumseq;
 611
 612                 calltmalign( *mseq1pt, *mseq2pt, *equivpt, (*strfiles)[i], (*chainids)[i], (*strfiles)[i], (*chainids)[i], *alloclenpt );
 613                 fprintf( stdout, ">%d_%s-%s\n%s\n", i+1, (*strfiles)[i], (*chainids)[i], (*seqpt)[i] );
 614                 alreadyoutput[i] = 1;
 615         }
 616 }
 617
 618 static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***seqpt, char ***mseq1pt, char ***mseq2pt, double **equivpt, int *alloclenpt )
 619 {
 620         int i, res;
 621         char *dumseq;
 622         char line[1000];
 623         char fname[1000];
 624         char command[1000];
 625         int linenum, istr, nstr;
 626         FILE *checkfp;
 627         char *sline;
 628         int use[1000];
 629         linenum = 0;
 630         nstr = 0;
 631         while( 1 )
 632         {
 633                 fgets( line, 999, fp );
 634                 if( feof( fp ) ) break;
 635                 sline = strip( line );
 636                 use[linenum] = 1;
 637                 if( sline[0] == '#' || strlen( sline ) < 2 )
 638                 {
 639                         use[linenum] = 0;
 640                         linenum++;
 641                         continue;
 642                 }
 643                 extractfirstword( sline );
 644                 checkfp = fopen( sline, "r" );
 645                 if( checkfp == NULL )
 646                 {
 647                         fprintf( stderr, "Cannot open %s.\n", sline );
 648                         exit( 1 );
 649                 }
 650 #if 0
 651                 fgets( linec, 999, checkfp );
 652                 if( strncmp( "HEADER ", linec, 7 ) )
 653                 {
 654                         fprintf( stderr, "Check the format of %s.\n", sline );
 655                         exit( 1 );
 656                 }
 657 #endif
 658                 if( checkcbeta( checkfp ) )
 659                 {
 660                         fprintf( stderr, "%s has no C-beta atoms.\n", sline );
 661                         exit( 1 );
 662                 }
 663                 else
 664                         nstr++;
 665                 fclose( checkfp );
 666                 linenum++;
 667         }
 668         njob = nstr;
 669         fprintf( stderr, "nstr = %d\n", nstr );
 670
 671         *strfiles = AllocateCharMtx( nstr, 1000 );
 672         *chainids = AllocateCharMtx( nstr, 2 );
 673
 674         rewind( fp );
 675         istr = 0;
 676         linenum = 0;
 677         while( 1 )
 678         {
 679                 fgets( line, 999, fp );
 680                 if( feof( fp ) ) break;
 681                 sline = strip( line );
 682                 if( use[linenum++] )
 683                 {
 684                         (*chainids)[istr][0] = getchainid( sline );
 685                         (*chainids)[istr][1] = 0;
 686                         extractfirstword( sline );
 687                         sprintf( fname, "%s", sline );
 688                         cutpath( fname );
 689                         sprintf( command, "cp %s %s.pdb", sline, fname );
 690                         system( command );
 691                         sprintf( command, "perl \"%s/clean.pl\" %s.pdb", whereispairalign, fname );
 692                         res = system( command );
 693                         if( res )
 694                         {
 695                                 fprintf( stderr, "error: Install clean.pl\n" );
 696                                 exit( 1 );
 697                         }
 698                         strcpy( (*strfiles)[istr++], fname );
 699                 }
 700         }
 701
 702         *seqpt = AllocateCharMtx( njob, nlenmax*2+1 );
 703         *mseq1pt = AllocateCharMtx( njob, 0 );
 704         *mseq2pt = AllocateCharMtx( njob, 0 );
 705         *equivpt = AllocateDoubleVec( nlenmax*2+1 );
 706         *alloclenpt = nlenmax*2;
 707         dumseq = AllocateCharVec( nlenmax*2+1 );
 708         alreadyoutput = AllocateIntVec( njob );
 709         for( i=0; i<njob; i++ ) alreadyoutput[i] = 0;
 710
 711         for( i=0; i<istr; i++ )
 712         {
 713                 fprintf( stderr, "i=%d\n", i );
 714                 (*seqpt)[i][0] = 0;
 715
 716                 (*mseq1pt)[0] = (*seqpt)[i];
 717                 (*mseq2pt)[0] = dumseq;
 718
 719                 callrash( i, i, *mseq1pt, *mseq2pt, *equivpt, (*strfiles)[i], (*chainids)[i], (*strfiles)[i], (*chainids)[i], *alloclenpt );
 720                 fprintf( stdout, ">%d_%s-%s\n%s\n", i+1, (*strfiles)[i], (*chainids)[i], (*seqpt)[i] );
 721                 alreadyoutput[i] = 1;
 722         }
 723 }
 724
 725 void arguments( int argc, char *argv[] )
 726 {
 727     int c;
 728
 729         usecache = 0;
 730         scale = 1.0;
 731         equivthreshold = 5;
 732         equivwinsize = 5;
 733         equivshortestlen = 1;
 734         inputfile = NULL;
 735         fftkeika = 0;
 736         pslocal = -1000.0;
 737         constraint = 0;
 738         nblosum = 62;
 739         fmodel = 0;
 740         calledByXced = 0;
 741         devide = 0;
 742         use_fft = 0;
 743         fftscore = 1;
 744         fftRepeatStop = 0;
 745         fftNoAnchStop = 0;
 746     weight = 3;
 747     utree = 1;
 748         tbutree = 1;
 749     refine = 0;
 750     check = 1;
 751     cut = 0.0;
 752     disp = 0;
 753     outgap = 1;
 754     alg = 'R';
 755     mix = 0;
 756         tbitr = 0;
 757         scmtd = 5;
 758         tbweight = 0;
 759         tbrweight = 3;
 760         checkC = 0;
 761         treemethod = 'x';
 762         contin = 0;
 763         scoremtx = 1;
 764         kobetsubunkatsu = 0;
 765         divpairscore = 0;
 766         dorp = NOTSPECIFIED;
 767         ppenalty = NOTSPECIFIED;
 768         ppenalty_OP = NOTSPECIFIED;
 769         ppenalty_ex = NOTSPECIFIED;
 770         ppenalty_EX = NOTSPECIFIED;
 771         poffset = NOTSPECIFIED;
 772         kimuraR = NOTSPECIFIED;
 773         pamN = NOTSPECIFIED;
 774         geta2 = GETA2;
 775         fftWinSize = NOTSPECIFIED;
 776         fftThreshold = NOTSPECIFIED;
 777         RNAppenalty = NOTSPECIFIED;
 778         RNApthr = NOTSPECIFIED;
 779
 780     while( --argc > 0 && (*++argv)[0] == '-' )
 781         {
 782         while ( ( c = *++argv[0] ) )
 783                 {
 784             switch( c )
 785             {
 786                                 case 'i':
 787                                         inputfile = *++argv;
 788                                         fprintf( stderr, "inputfile = %s\n", inputfile );
 789                                         --argc;
 790                                         goto nextoption;
 791                                 case 'f':
 792                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 793                                         --argc;
 794                                         goto nextoption;
 795                                 case 'g':
 796                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
 797                                         --argc;
 798                                         goto nextoption;
 799                                 case 'O':
 800                                         ppenalty_OP = (int)( atof( *++argv ) * 1000 - 0.5 );
 801                                         --argc;
 802                                         goto nextoption;
 803                                 case 'E':
 804                                         ppenalty_EX = (int)( atof( *++argv ) * 1000 - 0.5 );
 805                                         --argc;
 806                                         goto nextoption;
 807                                 case 'h':
 808                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
 809                                         --argc;
 810                                         goto nextoption;
 811                                 case 'k':
 812                                         kimuraR = atoi( *++argv );
 813 //                                      fprintf( stderr, "kimuraR = %d\n", kimuraR );
 814                                         --argc;
 815                                         goto nextoption;
 816                                 case 'b':
 817                                         nblosum = atoi( *++argv );
 818                                         scoremtx = 1;
 819 //                                      fprintf( stderr, "blosum %d\n", nblosum );
 820                                         --argc;
 821                                         goto nextoption;
 822                                 case 'j':
 823                                         pamN = atoi( *++argv );
 824                                         scoremtx = 0;
 825                                         TMorJTT = JTT;
 826                                         fprintf( stderr, "jtt %d\n", pamN );
 827                                         --argc;
 828                                         goto nextoption;
 829                                 case 'm':
 830                                         pamN = atoi( *++argv );
 831                                         scoremtx = 0;
 832                                         TMorJTT = TM;
 833                                         fprintf( stderr, "TM %d\n", pamN );
 834                                         --argc;
 835                                         goto nextoption;
 836                                 case 'd':
 837                                         whereispairalign = *++argv;
 838                                         fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
 839                                         --argc;
 840                                         goto nextoption;
 841                                 case 't':
 842                                         equivthreshold = atoi( *++argv );
 843                                         --argc;
 844                                         goto nextoption;
 845                                 case 'w':
 846                                         equivwinsize = atoi( *++argv );
 847                                         --argc;
 848                                         goto nextoption;
 849                                 case 'l':
 850                                         equivshortestlen = atoi( *++argv );
 851                                         --argc;
 852                                         goto nextoption;
 853                                 case 's':
 854                                         scale = atof( *++argv );
 855                                         --argc;
 856                                         goto nextoption;
 857                                 case 'c':
 858                                         usecache = 1;
 859                                         break;
 860 #if 1
 861                                 case 'a':
 862                                         fmodel = 1;
 863                                         break;
 864 #endif
 865                                 case 'r':
 866                                         fmodel = -1;
 867                                         break;
 868                                 case 'D':
 869                                         dorp = 'd';
 870                                         break;
 871                                 case 'P':
 872                                         dorp = 'p';
 873                                         break;
 874                                 case 'e':
 875                                         fftscore = 0;
 876                                         break;
 877 #if 0
 878                                 case 'O':
 879                                         fftNoAnchStop = 1;
 880                                         break;
 881 #endif
 882                                 case 'Q':
 883                                         calledByXced = 1;
 884                                         break;
 885                                 case 'x':
 886                                         disp = 1;
 887                                         break;
 888 #if 0
 889                                 case 'a':
 890                                         alg = 'a';
 891                                         break;
 892 #endif
 893                                 case 'S':
 894                                         alg = 'S';
 895                                         break;
 896                                 case 'L':
 897                                         alg = 'L';
 898                                         break;
 899                                 case 'B':
 900                                         alg = 'B';
 901                                         break;
 902                                 case 'T':
 903                                         alg = 'T';
 904                                         break;
 905                                 case 'H':
 906                                         alg = 'H';
 907                                         break;
 908                                 case 'M':
 909                                         alg = 'M';
 910                                         break;
 911                                 case 'R':
 912                                         alg = 'R';
 913                                         break;
 914                                 case 'N':
 915                                         alg = 'N';
 916                                         break;
 917                                 case 'K':
 918                                         alg = 'K';
 919                                         break;
 920                                 case 'A':
 921                                         alg = 'A';
 922                                         break;
 923                                 case 'V':
 924                                         alg = 'V';
 925                                         break;
 926                                 case 'C':
 927                                         alg = 'C';
 928                                         break;
 929                                 case 'F':
 930                                         use_fft = 1;
 931                                         break;
 932                                 case 'v':
 933                                         tbrweight = 3;
 934                                         break;
 935                                 case 'y':
 936                                         divpairscore = 1;
 937                                         break;
 938 /* Modified 01/08/27, default: user tree */
 939                                 case 'J':
 940                                         tbutree = 0;
 941                                         break;
 942 /* modification end. */
 943 #if 0
 944                                 case 'z':
 945                                         fftThreshold = atoi( *++argv );
 946                                         --argc;
 947                                         goto nextoption;
 948                                 case 'w':
 949                                         fftWinSize = atoi( *++argv );
 950                                         --argc;
 951                                         goto nextoption;
 952                                 case 'Z':
 953                                         checkC = 1;
 954                                         break;
 955 #endif
 956                 default:
 957                     fprintf( stderr, "illegal option %c\n", c );
 958                     argc = 0;
 959                     break;
 960             }
 961                 }
 962                 nextoption:
 963                         ;
 964         }
 965     if( argc == 1 )
 966     {
 967         cut = atof( (*argv) );
 968         argc--;
 969     }
 970     if( argc != 0 )
 971     {
 972         fprintf( stderr, "options: Check source file !\n" );
 973         exit( 1 );
 974     }
 975         if( tbitr == 1 && outgap == 0 )
 976         {
 977                 fprintf( stderr, "conflicting options : o, m or u\n" );
 978                 exit( 1 );
 979         }
 980         if( alg == 'C' && outgap == 0 )
 981         {
 982                 fprintf( stderr, "conflicting options : C, o\n" );
 983                 exit( 1 );
 984         }
 985 }
 986
 987 int countamino( char *s, int end )
 988 {
 989         int val = 0;
 990         while( end-- )
 991                 if( *s++ != '-' ) val++;
 992         return( val );
 993 }
 994
 995 static void pairalign( char name[M][B], int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, double *equiv, double *effarr, char **strfiles, char **chainids, int alloclen )
 996 {
 997         int i, j, ilim;
 998         int clus1, clus2;
 999         int off1, off2;
1000         float pscore = 0.0; // by D.Mathog
1001         static char *indication1, *indication2;
1002         FILE *hat2p, *hat3p;
1003         static double **distancemtx;
1004         static double *effarr1 = NULL;
1005         static double *effarr2 = NULL;
1006         char *pt;
1007         char *hat2file = "hat2";
1008         LocalHom **localhomtable, *tmpptr;
1009         static char **pair;
1010 //      int intdum;
1011         double bunbo;
1012         char **checkseq;
1013
1014
1015         localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
1016         for( i=0; i<njob; i++)
1017         {
1018
1019                 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1020                 for( j=0; j<njob; j++)
1021                 {
1022                         localhomtable[i][j].start1 = -1;
1023                         localhomtable[i][j].end1 = -1;
1024                         localhomtable[i][j].start2 = -1;
1025                         localhomtable[i][j].end2 = -1;
1026                         localhomtable[i][j].opt = -1.0;
1027                         localhomtable[i][j].next = NULL;
1028                         localhomtable[i][j].nokori = 0;
1029                 }
1030         }
1031
1032         if( effarr1 == NULL )
1033         {
1034                 distancemtx = AllocateDoubleMtx( njob, njob );
1035                 effarr1 = AllocateDoubleVec( njob );
1036                 effarr2 = AllocateDoubleVec( njob );
1037                 indication1 = AllocateCharVec( 150 );
1038                 indication2 = AllocateCharVec( 150 );
1039                 checkseq = AllocateCharMtx( njob, alloclen );
1040 #if 0
1041 #else
1042                 pair = AllocateCharMtx( njob, njob );
1043 #endif
1044         }
1045
1046 #if 0
1047         fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
1048 #endif
1049
1050 #if 0
1051         for( i=0; i<njob; i++ )
1052                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
1053 #endif
1054
1055
1056 //      writePre( njob, name, nlen, aseq, 0 );
1057
1058         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
1059         for( i=0; i<njob; i++ ) pair[i][i] = 1;
1060
1061         for( i=0; i<njob; i++ )
1062         {
1063                 strcpy( checkseq[i], seq[i] );
1064 //              fprintf( stderr, "checkseq[%d] = %s\n", i, checkseq[i] );
1065         }
1066
1067
1068         ilim = njob - 1;
1069         for( i=0; i<ilim; i++ )
1070         {
1071                 fprintf( stderr, "% 5d / %d\r", i, njob );
1072
1073
1074                 for( j=i+1; j<njob; j++ )
1075                 {
1076
1077
1078 #if 0
1079                         if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
1080                         {
1081                                 distancemtx[i][j] = pscore;
1082                                 continue;
1083                         }
1084 #endif
1085
1086                         strcpy( aseq[i], seq[i] );
1087                         strcpy( aseq[j], seq[j] );
1088                         clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
1089                         clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
1090         //              fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
1091         //              fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
1092
1093 #if 0
1094                         fprintf( stderr, "group1 = %.66s", indication1 );
1095                         fprintf( stderr, "\n" );
1096                         fprintf( stderr, "group2 = %.66s", indication2 );
1097                         fprintf( stderr, "\n" );
1098 #endif
1099 //                      for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
1100
1101 #if 1
1102                         {
1103                                 switch( alg )
1104                                 {
1105                                         case( 'T' ):
1106                                                 fprintf( stderr, "  Calling tmalign %d-%d/%d    \r", i+1, j+1, njob );
1107                                                 pscore = calltmalign( mseq1, mseq2, equiv, strfiles[i], chainids[i], strfiles[j], chainids[j], alloclen );
1108                                                 off1 = off2 = 0;
1109                                                 break;
1110                                         case( 'R' ):
1111                                                 fprintf( stderr, "  Calling PDP_ASH.pl %d-%d/%d    \r", i+1, j+1, njob );
1112                                                 pscore = callrash( i, j, mseq1, mseq2, equiv, strfiles[i], chainids[i], strfiles[j], chainids[j], alloclen );
1113                                                 off1 = off2 = 0;
1114                                                 break;
1115                                         ErrorExit( "ERROR IN SOURCE FILE" );
1116                                 }
1117                         }
1118 #endif
1119                         distancemtx[i][j] = pscore;
1120 #if SCOREOUT
1121                         fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
1122 #endif
1123
1124                         putlocalhom_str( mseq1[0], mseq2[0], equiv, scale, localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
1125 #if 1
1126
1127                         if( alreadyoutput[i] == 0 )
1128                         {
1129                                 alreadyoutput[i] = 1;
1130                                 gappick0( seq[i], mseq1[0] );
1131                                 fprintf( stdout, ">%d_%s-%s\n%s\n", i+1, strfiles[i], chainids[i], seq[i] );
1132                                 strcpy( checkseq[i], seq[i] );
1133                         }
1134                         else
1135                         {
1136                                 gappick0( seq[i], mseq1[0] );
1137                                 fprintf( stderr, "checking seq%d\n", i );
1138
1139 //                              fprintf( stderr, "     seq=%s\n", seq[i] );
1140 //                              fprintf( stderr, "checkseq=%s\n", checkseq[i] );
1141
1142                                 if( strcmp( checkseq[i], seq[i] ) )
1143                                 {
1144                                         fprintf( stderr, "\n\nWARNING: Sequence changed!!\n" );
1145                                         fprintf( stderr, "i=%d\n", i );
1146                                         fprintf( stderr, "     seq=%s\n", seq[i] );
1147                                         fprintf( stderr, "checkseq=%s\n", checkseq[i] );
1148                                         exit( 1 );
1149                                 }
1150                         }
1151                         if( alreadyoutput[j] == 0 )
1152                         {
1153                                 alreadyoutput[j] = 1;
1154                                 gappick0( seq[j], mseq2[0] );
1155                                 fprintf( stdout, ">%d_%s-%s\n%s\n", j+1, strfiles[j], chainids[j], seq[j] );
1156                                 strcpy( checkseq[j], seq[j] );
1157                         }
1158                         else
1159                         {
1160                                 gappick0( seq[j], mseq2[0] );
1161                                 fprintf( stderr, "checking seq%d\n", j );
1162                                 if( strcmp( checkseq[j], seq[j] ) )
1163                                 {
1164                                         fprintf( stderr, "\n\nWARNING: Sequence changed!!\n" );
1165                                         fprintf( stderr, "j=%d\n", j );
1166                                         fprintf( stderr, "     seq=%s\n", seq[j] );
1167                                         fprintf( stderr, "checkseq=%s\n", checkseq[j] );
1168                                         exit( 1 );
1169                                 }
1170                         }
1171 #endif
1172                 }
1173         }
1174         for( i=0; i<njob; i++ )
1175         {
1176                 pscore = 0.0;
1177                 for( pt=seq[i]; *pt; pt++ )
1178                         pscore += amino_dis[(int)*pt][(int)*pt];
1179                 distancemtx[i][i] = pscore;
1180
1181         }
1182
1183         ilim = njob-1;
1184         for( i=0; i<ilim; i++ )
1185         {
1186                 for( j=i+1; j<njob; j++ )
1187                 {
1188                         bunbo = MIN( distancemtx[i][i], distancemtx[j][j] );
1189                         if( bunbo == 0.0 )
1190                                 distancemtx[i][j] = 2.0;
1191                         else
1192                                 distancemtx[i][j] = ( 1.0 - distancemtx[i][j] / bunbo ) * 2.0;
1193                 }
1194         }
1195
1196         hat2p = fopen( hat2file, "w" );
1197         if( !hat2p ) ErrorExit( "Cannot open hat2." );
1198         WriteHat2( hat2p, njob, name, distancemtx );
1199         fclose( hat2p );
1200
1201         fprintf( stderr, "##### writing hat3\n" );
1202         hat3p = fopen( "hat3", "w" );
1203         if( !hat3p ) ErrorExit( "Cannot open hat3." );
1204         ilim = njob-1;
1205         for( i=0; i<ilim; i++ )
1206         {
1207                 for( j=i+1; j<njob; j++ )
1208                 {
1209                         for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
1210                         {
1211                                 if( tmpptr->opt == -1.0 ) continue;
1212                                 fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d k\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2 );
1213                         }
1214                 }
1215         }
1216         fclose( hat3p );
1217 #if DEBUG
1218         fprintf( stderr, "calling FreeLocalHomTable\n" );
1219 #endif
1220         FreeLocalHomTable( localhomtable, njob );
1221 #if DEBUG
1222         fprintf( stderr, "done. FreeLocalHomTable\n" );
1223 #endif
1224 }
1225
1226 static void WriteOptions( FILE *fp )
1227 {
1228
1229         if( dorp == 'd' ) fprintf( fp, "DNA\n" );
1230         else
1231         {
1232                 if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
1233                 else if( scoremtx ==  1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
1234                 else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
1235         }
1236     fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1237     if( use_fft ) fprintf( fp, "FFT on\n" );
1238
1239         fprintf( fp, "tree-base method\n" );
1240         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
1241         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
1242         if( tbitr || tbweight )
1243         {
1244                 fprintf( fp, "iterate at each step\n" );
1245                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
1246                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
1247                 if( tbweight ) fprintf( fp, "  weighted\n" );
1248                 fprintf( fp, "\n" );
1249         }
1250
1251          fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1252
1253         if( alg == 'a' )
1254                 fprintf( fp, "Algorithm A\n" );
1255         else if( alg == 'A' )
1256                 fprintf( fp, "Algorithm A+\n" );
1257         else if( alg == 'S' )
1258                 fprintf( fp, "Apgorithm S\n" );
1259         else if( alg == 'C' )
1260                 fprintf( fp, "Apgorithm A+/C\n" );
1261         else
1262                 fprintf( fp, "Unknown algorithm\n" );
1263
1264     if( use_fft )
1265     {
1266         fprintf( fp, "FFT on\n" );
1267         if( dorp == 'd' )
1268             fprintf( fp, "Basis : 4 nucleotides\n" );
1269         else
1270         {
1271             if( fftscore )
1272                 fprintf( fp, "Basis : Polarity and Volume\n" );
1273             else
1274                 fprintf( fp, "Basis : 20 amino acids\n" );
1275         }
1276         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
1277         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
1278     }
1279         else
1280         fprintf( fp, "FFT off\n" );
1281         fflush( fp );
1282 }
1283
1284
1285 int main( int argc, char *argv[] )
1286 {
1287         static int  nlen[M];
1288         static char name[M][B], **seq;
1289         static char **mseq1, **mseq2;
1290         static char **aseq;
1291         static char **bseq;
1292         static double *eff;
1293         static double *equiv;
1294         char **strfiles;
1295         char **chainids;
1296         int i;
1297         FILE *infp;
1298         char c;
1299         int alloclen;
1300
1301         arguments( argc, argv );
1302
1303         if( equivthreshold < 1 || 9 < equivthreshold )
1304         {
1305                 fprintf( stderr, "-t n, n must be 1..9\n" );
1306                 exit( 1 );
1307         }
1308
1309         if( ( equivwinsize + 1 ) % 2 != 0 )
1310         {
1311                 fprintf( stderr, "equivwinsize = %d\n", equivwinsize );
1312                 fprintf( stderr, "It must be an odd number.\n" );
1313                 exit( 1 );
1314         }
1315
1316         if( inputfile )
1317         {
1318                 infp = fopen( inputfile, "r" );
1319                 if( !infp )
1320                 {
1321                         fprintf( stderr, "Cannot open %s\n", inputfile );
1322                         exit( 1 );
1323                 }
1324         }
1325         else
1326                 infp = stdin;
1327
1328         nlenmax = 10000; // tekitou
1329
1330         if( alg == 'R' )
1331                 prepareash( infp, &strfiles, &chainids, &seq, &mseq1, &mseq2, &equiv, &alloclen );
1332         else if( alg == 'T' )
1333                 preparetmalign( infp, &strfiles, &chainids, &seq, &mseq1, &mseq2, &equiv, &alloclen );
1334
1335         fclose( infp );
1336
1337         aseq = AllocateCharMtx( njob, nlenmax*2+1 );
1338         bseq = AllocateCharMtx( njob, nlenmax*2+1 );
1339         eff = AllocateDoubleVec( njob );
1340
1341         for( i=0; i<njob; i++ )
1342         {
1343                 fprintf( stderr, "str%d = %s-%s\n", i, strfiles[i], chainids[i] );
1344         }
1345
1346         if( njob < 1 )
1347         {
1348                 fprintf( stderr, "No structure found.\n" );
1349                 exit( 1 );
1350         }
1351         if( njob < 2 )
1352         {
1353                 fprintf( stderr, "Only %d structure found.\n", njob );
1354                 exit( 0 );
1355         }
1356         if( njob > M )
1357         {
1358                 fprintf( stderr, "The number of structures must be < %d\n", M );
1359                 fprintf( stderr, "Please try sequence-based methods for such large data.\n" );
1360                 exit( 1 );
1361         }
1362
1363
1364
1365 #if 0
1366         readData( infp, name, nlen, seq );
1367 #endif
1368
1369         constants( njob, seq );
1370
1371 #if 0
1372         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
1373 #endif
1374
1375         initSignalSM();
1376
1377         initFiles();
1378
1379         WriteOptions( trap_g );
1380
1381         c = seqcheck( seq );
1382         if( c )
1383         {
1384                 fprintf( stderr, "Illegal character %c\n", c );
1385                 exit( 1 );
1386         }
1387
1388 //      writePre( njob, name, nlen, seq, 0 );
1389
1390         for( i=0; i<njob; i++ ) eff[i] = 1.0;
1391
1392
1393         for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1394
1395         pairalign( name, nlen, bseq, aseq, mseq1, mseq2, equiv, eff, strfiles, chainids, alloclen );
1396
1397         fprintf( trap_g, "done.\n" );
1398 #if DEBUG
1399         fprintf( stderr, "closing trap_g\n" );
1400 #endif
1401         fclose( trap_g );
1402
1403 //      writePre( njob, name, nlen, aseq, !contin );
1404 #if 0
1405         writeData( stdout, njob, name, nlen, aseq );
1406 #endif
1407 #if IODEBUG
1408         fprintf( stderr, "OSHIMAI\n" );
1409 #endif
1410         SHOWVERSION;
1411         return( 0 );
1412 }