binaries/src/mafft/core/pairlocalalign.c

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