binaries/src/mafft/core/dvtditr.c

   1  /* Tree-dependent-iteration */
   2  /* Devide to segments       */
   3
   4 #include "mltaln.h"
   5
   6 extern char **seq_g;
   7 extern char **res_g;
   8 static int subalignment;
   9 static int subalignmentoffset;
  10 static int specifictarget;
  11
  12 static int intop;
  13 static int intree;
  14 static double autosubalignment;
  15
  16
  17 static void calcmaxdistclass( void )
  18 {
  19         int c;
  20         double rep;
  21         for( c=0; c<ndistclass; c++ )
  22         {
  23                 rep = (double) 2 * c / ndistclass; // dist:0-2 for dist2offset
  24 //              fprintf( stderr, "c=%d, rep=%f, offset=%f\n", c, rep, dist2offset( rep )  );
  25                 if( dist2offset( rep ) == 0.0 )
  26                         break;
  27         }
  28         fprintf( stderr, "ndistclass = %d, maxdistclass = %d\n", ndistclass, c+1 );
  29         maxdistclass = c + 1;
  30 //      maxdistclass = ndistclass; // CHUUI!!!!
  31         return;
  32 }
  33
  34 void arguments( int argc, char *argv[] )
  35 {
  36         int c;
  37         char *argkey;
  38
  39         outnumber = 0;
  40         nthread = 1;
  41         randomseed = 0;
  42         scoreout = 0;
  43         spscoreout = 0;
  44         parallelizationstrategy = BAATARI1;
  45         intop = 0;
  46         intree = 0;
  47         inputfile = NULL;
  48         rnakozo = 0;
  49         rnaprediction = 'm';
  50         nevermemsave = 0;
  51         score_check = 1;
  52         fftkeika = 1;
  53         constraint = 0;
  54         fmodel = 0;
  55         kobetsubunkatsu = 1;
  56         bunkatsu = 1;
  57         nblosum = 62;
  58         niter = 100;
  59         calledByXced = 0;
  60         devide = 1;
  61         divWinSize = 20; /* 70 */
  62         divThreshold = 65;
  63         fftscore = 1;
  64         fftRepeatStop = 0;
  65         fftNoAnchStop = 0;
  66     scmtd = 5;
  67         cooling = 1;
  68     weight = 4;
  69     utree = 1;
  70     refine = 1;
  71     check = 1;
  72     cut = 0.0;
  73         disp = 0;
  74         outgap = 1;
  75         use_fft = 0; // CHUUI dochira demo mafft.tmpl deha F
  76         force_fft = 0;
  77         alg = 'A';  /* chuui */
  78         mix = 0;
  79         checkC = 0;
  80         tbitr = 0;
  81         treemethod = 'X';
  82         sueff_global = 0.1;
  83         scoremtx = 1;
  84         dorp = NOTSPECIFIED;
  85         ppenalty = NOTSPECIFIED;
  86         penalty_shift_factor = 1000.0;
  87         ppenalty_ex = NOTSPECIFIED;
  88         poffset = NOTSPECIFIED;
  89         kimuraR = NOTSPECIFIED;
  90         pamN = NOTSPECIFIED;
  91         geta2 = GETA2;
  92         fftWinSize = NOTSPECIFIED;
  93         fftThreshold = NOTSPECIFIED;
  94         RNAppenalty = NOTSPECIFIED;
  95         RNAppenalty_ex = NOTSPECIFIED;
  96         RNApthr = NOTSPECIFIED;
  97         TMorJTT = JTT;
  98         consweight_multi = 1.0;
  99         consweight_rna = 0.0;
 100         subalignment = 0;
 101         subalignmentoffset = 0;
 102         legacygapcost = 0;
 103         specificityconsideration = 0.0;
 104         autosubalignment = 0.0;
 105         specifictarget = 0;
 106         nwildcard = 0;
 107
 108         while( --argc > 0 && (*++argv)[0] == '-' )
 109         {
 110                 while ( (c = *++argv[0]) )
 111                 {
 112                         switch( c )
 113                         {
 114                                 case 'i':
 115                                         inputfile = *++argv;
 116                                         fprintf( stderr, "inputfile = %s\n", inputfile );
 117                                         --argc;
 118                                         goto nextoption;
 119                                 case 'I':
 120                                         niter = myatoi( *++argv );
 121                                         fprintf( stderr, "niter = %d\n", niter );
 122                                         --argc;
 123                                         goto nextoption;
 124                                 case 'e':
 125                                         RNApthr = (int)( atof( *++argv ) * 1000 - 0.5 );
 126                                         --argc;
 127                                         goto nextoption;
 128                                 case 'o':
 129                                         RNAppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 130                                         --argc;
 131                                         goto nextoption;
 132                                 case 'f':
 133                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 134 //                                      fprintf( stderr, "ppenalty = %d\n", ppenalty );
 135                                         --argc;
 136                                         goto nextoption;
 137                                 case 'Q':
 138                                         penalty_shift_factor = atof( *++argv );
 139                                         if( penalty_shift_factor < 100.0 && penalty_shift_factor != 2.0 )
 140                                         {
 141                                                 fprintf( stderr, "%f, penalty_shift is fixed to penalty x 2 in the iterative refinement phase.\n", penalty_shift_factor );
 142                                                 penalty_shift_factor = 2.0;
 143                                         }
 144                                         --argc;
 145                                         goto nextoption;
 146                                 case 'g':
 147                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
 148 //                                      fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
 149                                         --argc;
 150                                         goto nextoption;
 151                                 case 'h':
 152                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
 153                                         fprintf( stderr, "poffset = %d\n", poffset );
 154                                         --argc;
 155                                         goto nextoption;
 156                                 case 'k':
 157                                         kimuraR = myatoi( *++argv );
 158                                         fprintf( stderr, "kappa = %d\n", kimuraR );
 159                                         --argc;
 160                                         goto nextoption;
 161                                 case 'b':
 162                                         nblosum = myatoi( *++argv );
 163                                         scoremtx = 1;
 164                                         fprintf( stderr, "blosum %d / kimura 200\n", nblosum );
 165                                         --argc;
 166                                         goto nextoption;
 167                                 case 'j':
 168                                         pamN = myatoi( *++argv );
 169                                         scoremtx = 0;
 170                                         TMorJTT = JTT;
 171                                         fprintf( stderr, "jtt/kimura %d\n", pamN );
 172                                         --argc;
 173                                         goto nextoption;
 174                                 case 'm':
 175                                         pamN = myatoi( *++argv );
 176                                         scoremtx = 0;
 177                                         TMorJTT = TM;
 178                                         fprintf( stderr, "tm %d\n", pamN );
 179                                         --argc;
 180                                         goto nextoption;
 181                                 case 'l':
 182                                         fastathreshold = atof( *++argv );
 183                                         constraint = 2;
 184                                         --argc;
 185                                         goto nextoption;
 186                                 case 'r':
 187                                         consweight_rna = atof( *++argv );
 188                                         rnakozo = 1;
 189                                         --argc;
 190                                         goto nextoption;
 191                                 case 'c':
 192                                         consweight_multi = atof( *++argv );
 193                                         --argc;
 194                                         goto nextoption;
 195                                 case 'C':
 196                                         nthread = myatoi( *++argv );
 197                                         fprintf( stderr, "nthread = %d\n", nthread );
 198                                         --argc;
 199                                         goto nextoption;
 200                                 case 'H':
 201                                         subalignment = 1;
 202                                         subalignmentoffset = myatoi( *++argv );
 203                                         --argc;
 204                                         goto nextoption;
 205                                 case 't':
 206                                         randomseed = myatoi( *++argv );
 207                                         fprintf( stderr, "randomseed = %d\n", randomseed );
 208                                         --argc;
 209                                         goto nextoption;
 210                                 case 'p':
 211                                         argkey = *++argv;
 212                                         if( !strcmp( argkey, "BESTFIRST" ) ) parallelizationstrategy = BESTFIRST;
 213                                         else if( !strcmp( argkey, "BAATARI0" ) ) parallelizationstrategy = BAATARI0;
 214                                         else if( !strcmp( argkey, "BAATARI1" ) ) parallelizationstrategy = BAATARI1;
 215                                         else if( !strcmp( argkey, "BAATARI2" ) ) parallelizationstrategy = BAATARI2;
 216                                         else
 217                                         {
 218                                                 fprintf( stderr, "Unknown parallelization strategy, %s\n", argkey );
 219                                                 exit( 1 );
 220                                         }
 221 //                                      exit( 1 );
 222                                         --argc;
 223                                         goto nextoption;
 224                                 case 's':
 225                                         specificityconsideration = (double)myatof( *++argv );
 226 //                                      fprintf( stderr, "specificityconsideration = %f\n", specificityconsideration );
 227                                         --argc;
 228                                         goto nextoption;
 229 #if 0
 230                                 case 'S' :
 231                                         scoreout = 1; // for checking parallel calculation
 232                                         break;
 233 #else
 234                                 case 'S' :
 235                                         spscoreout = 1; // 2014/Dec/30, sp score
 236                                         break;
 237 #endif
 238 #if 0
 239                                 case 's' :
 240                                         RNAscoremtx = 'r';
 241                                         break;
 242 #endif
 243 #if 1
 244                                 case 'a':
 245                                         fmodel = 1;
 246                                         break;
 247 #endif
 248                                 case 'N':
 249                                         nevermemsave = 1;
 250                                         break;
 251                                 case 'D':
 252                                         dorp = 'd';
 253                                         break;
 254                                 case 'P':
 255                                         dorp = 'p';
 256                                         break;
 257 #if 0
 258                                 case 'Q':
 259                                         alg = 'Q';
 260                                         break;
 261 #endif
 262                                 case 'R':
 263                                         rnaprediction = 'r';
 264                                         break;
 265                                 case 'O':
 266                                         fftNoAnchStop = 1;
 267                                         break;
 268 #if 0
 269                                 case 'e':
 270                                         fftscore = 0;
 271                                         break;
 272                                 case 'r':
 273                                         fmodel = -1;
 274                                         break;
 275                                 case 'R':
 276                                         fftRepeatStop = 1;
 277                                         break;
 278 #endif
 279                                 case 'T':
 280                                         kobetsubunkatsu = 0;
 281                                         break;
 282                                 case 'B':
 283                                         bunkatsu = 0;
 284                                         break;
 285 #if 0
 286                                 case 'c':
 287                                         cooling = 1;
 288                                         break;
 289                                 case 'a':
 290                                         alg = 'a';
 291                                         break;
 292                                 case 's' :
 293                                         treemethod = 's';
 294                                         break;
 295                                 case 'H':
 296                                         alg = 'H';
 297                                         break;
 298 #endif
 299                                 case 'A':
 300                                         alg = 'A';
 301                                         break;
 302                                 case 'M':
 303                                         alg = 'M';
 304                                         break;
 305                                 case '@':
 306                                         alg = 'd';
 307                                         break;
 308                                 case 'F':
 309                                         use_fft = 1;
 310                                         break;
 311 #if 0
 312                                 case 't':
 313                                         weight = 4;
 314                                         break;
 315 #endif
 316                                 case 'u':
 317                                         weight = 0;
 318                                         break;
 319                                 case 'U':
 320                                         intree = 1;
 321                                         break;
 322                                 case 'V':
 323                                         intop = 1;
 324                                         break;
 325                                 case 'J':
 326                                         utree = 0;
 327                                         break;
 328 #if 0
 329                                 case 'd':
 330                                         disp = 1;
 331                                         break;
 332 #endif
 333                                 case 'Z':
 334                                         score_check = 0;
 335                                         break;
 336                                 case 'Y':
 337                                         score_check = 2;
 338                                         break;
 339                                 case 'L':
 340                                         legacygapcost = 1;
 341                                         break;
 342 #if 0
 343                                 case 'n' :
 344                                         treemethod = 'n';
 345                                         break;
 346 #endif
 347                                 case 'n' :
 348                                         outnumber = 1;
 349                                         break;
 350                                 case 'X':
 351                                         treemethod = 'X';
 352                                         sueff_global = atof( *++argv );
 353                                         fprintf( stderr, "sueff_global = %f\n", sueff_global );
 354                                         --argc;
 355                                         goto nextoption;
 356 #if 0
 357                                 case 'E' :
 358                                         treemethod = 'E';
 359                                         break;
 360                                 case 'q' :
 361                                         treemethod = 'q';
 362                                         break;
 363 #endif
 364                                 case 'E':
 365                                         autosubalignment = atof( *++argv );
 366                                         fprintf( stderr, "autosubalignment = %f\n", autosubalignment );
 367                                         --argc;
 368                                         goto nextoption;
 369                                 case 'W':
 370                                         minimumweight = atof( *++argv );
 371                                         fprintf( stderr, "minimumweight = %f\n", minimumweight );
 372                                         --argc;
 373                                         goto nextoption;
 374                                 case 'z':
 375                                         fftThreshold = myatoi( *++argv );
 376                                         --argc;
 377                                         goto nextoption;
 378                                 case 'w':
 379                                         fftWinSize = myatoi( *++argv );
 380                                         --argc;
 381                                         goto nextoption;
 382                                 case '=':
 383                                         specifictarget = 1;
 384                                         break;
 385                                 case ':':
 386                                         nwildcard = 1;
 387                                         break;
 388                                 default:
 389                                         fprintf( stderr, "illegal option %c\n", c );
 390                                         argc = 0;
 391                                         break;
 392                         }
 393                 }
 394                 nextoption:
 395                         ;
 396         }
 397         if( argc == 1 )
 398         {
 399                 cut = atof( (*argv) );
 400                 argc--;
 401         }
 402         if( argc != 0 )
 403         {
 404                 fprintf( stderr, "options : Check source file!\n" );
 405                 exit( 1 );
 406         }
 407 #if 0
 408         if( alg == 'A' && weight == 0 )
 409                 ErrorExit( "ERROR : Algorithm A+ and un-weighted\n" );
 410 #endif
 411 }
 412
 413
 414 int main( int argc, char *argv[] )
 415 {
 416     int identity;
 417         static int nlen[M];
 418         static char **name, **seq, **aseq, **bseq;
 419         static Segment *segment = NULL;
 420         static int anchors[MAXSEG];
 421         int i, j;
 422         int iseg, nseg;
 423         int ***topol;
 424         double **len;
 425         double **eff;
 426         FILE *prep;
 427         FILE *infp;
 428         FILE *orderfp;
 429         int alloclen;
 430         int returnvalue;
 431         char c;
 432         int ocut;
 433         char **seq_g_bk;
 434         LocalHom **localhomtable = NULL; // by D.Mathog
 435         RNApair ***singlerna;
 436         int nogaplen;
 437         static char **nogap1seq;
 438         static char *kozoarivec;
 439         int nkozo;
 440         int alignmentlength;
 441         int **skipthisbranch;
 442         int foundthebranch;
 443         int nsubalignments, maxmem;
 444         int **subtable;
 445         int *insubtable;
 446         int *preservegaps;
 447         char ***subalnpt;
 448         int ntarget, *targetmap, *targetmapr;
 449         int ilim;
 450
 451         arguments( argc, argv );
 452 #ifndef enablemultithread
 453         nthread = 0;
 454 #endif
 455         if( fastathreshold < 0.0001 ) constraint = 0;
 456
 457         if( inputfile )
 458         {
 459                 infp = fopen( inputfile, "r" );
 460                 if( !infp )
 461                 {
 462                         fprintf( stderr, "Cannot open %s\n", inputfile );
 463                         exit( 1 );
 464                 }
 465         }
 466         else
 467                 infp = stdin;
 468
 469 #if 0
 470         PreRead( stdin, &njob, &nlenmax );
 471 #else
 472         getnumlen( infp );
 473 #endif
 474         rewind( infp );
 475
 476         nkozo = 0;
 477
 478         if( njob < 2 )
 479         {
 480                 fprintf( stderr, "At least 2 sequences should be input!\n"
 481                                                  "Only %d sequence found.\n", njob );
 482                 exit( 1 );
 483         }
 484
 485
 486         if( subalignment )
 487         {
 488                 readsubalignmentstable( njob, NULL, NULL, &nsubalignments, &maxmem );
 489                 fprintf( stderr, "nsubalignments = %d\n", nsubalignments );
 490                 fprintf( stderr, "maxmem = %d\n", maxmem );
 491                 subtable = AllocateIntMtx( nsubalignments, maxmem+1 );
 492                 insubtable = AllocateIntVec( njob );
 493                 preservegaps = AllocateIntVec( njob );
 494                 for( i=0; i<njob; i++ ) insubtable[i] = 0;
 495                 for( i=0; i<njob; i++ ) preservegaps[i] = 0;
 496                 subalnpt = AllocateCharCub( nsubalignments, maxmem, 0 );
 497                 readsubalignmentstable( njob, subtable, preservegaps, NULL, NULL );
 498                 for( i=0; i<nsubalignments; i++ ) for( j=0; j<insubtable[i]; j++ )
 499                 {
 500                         if( subtable[i][j] < 0 )
 501                         {
 502                                 fprintf( stderr, "Not supported in the iterative refinmenment mode.\n" );
 503                                 fprintf( stderr, "Please use a positive number to represent a sequence.\n" );
 504                         }
 505                 }
 506         }
 507
 508         ocut = cut;
 509
 510         segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 511 //      if( treemethod == 'X' || treemethod == 'E' || treemethod == 'q' )
 512         topol = AllocateIntCub( njob, 2, njob );
 513         len = AllocateDoubleMtx( njob, 2 );
 514         eff = AllocateDoubleMtx( njob, njob );
 515         seq = AllocateCharMtx( njob, nlenmax*9+1 );
 516         name = AllocateCharMtx( njob, B+1 );
 517         seq_g = AllocateCharMtx( njob, nlenmax*9+1 );
 518         res_g = AllocateCharMtx( njob, nlenmax*9+1 );
 519         aseq = AllocateCharMtx( njob, nlenmax*9+1 );
 520         bseq = AllocateCharMtx( njob, nlenmax*9+1 );
 521         nogap1seq = AllocateCharMtx( 1, nlenmax*1+1 );
 522         alloclen = nlenmax * 9;
 523         seq_g_bk = AllocateCharMtx( njob, 0 );
 524         for( i=0; i<njob; i++ ) seq_g_bk[i] = seq_g[i];
 525         kozoarivec = AllocateCharVec( njob );
 526         skipthisbranch = AllocateIntMtx( njob, 2 );
 527         for( i=0; i<njob; i++ ) skipthisbranch[i][0] = skipthisbranch[i][1] = 0;
 528
 529
 530 #if 0
 531         Read( name, nlen, seq_g );
 532 #else
 533         readData_pointer( infp, name, nlen, seq_g );
 534 #endif
 535         fclose( infp );
 536
 537         if( specifictarget )
 538         {
 539                 targetmap = calloc( njob, sizeof( int ) );
 540                 ntarget = 0;
 541                 for( i=0; i<njob; i++ )
 542                 {
 543                         targetmap[i] = -1;
 544                         if( !strncmp( name[i]+1, "_focus_", 7 ) )
 545                                 targetmap[i] = ntarget++;
 546                 }
 547                 targetmapr = calloc( ntarget, sizeof( int ) );
 548                 for( i=0; i<njob; i++ )
 549                         if( targetmap[i] != -1 ) targetmapr[targetmap[i]] = i;
 550         }
 551         else
 552         {
 553                 ntarget = njob;
 554                 targetmap = calloc( njob, sizeof( int ) );
 555                 targetmapr = calloc( njob, sizeof( int ) );
 556                 for( i=0; i<njob; i++ )
 557                         targetmap[i] = targetmapr[i] = i;
 558         }
 559
 560         if( constraint )
 561         {
 562                 ilim = njob;
 563                 localhomtable = (LocalHom **)calloc( ntarget, sizeof( LocalHom *) );
 564                 for( i=0; i<ntarget; i++)
 565                 {
 566                         localhomtable[i] = (LocalHom *)calloc( ilim, sizeof( LocalHom ) );
 567                         for( j=0; j<ilim; j++)
 568                         {
 569                                 localhomtable[i][j].start1 = -1;
 570                                 localhomtable[i][j].end1 = -1;
 571                                 localhomtable[i][j].start2 = -1;
 572                                 localhomtable[i][j].end2 = -1;
 573                                 localhomtable[i][j].overlapaa = -1.0;
 574                                 localhomtable[i][j].opt = -1.0;
 575                                 localhomtable[i][j].importance = -1.0;
 576                                 localhomtable[i][j].next = NULL;
 577                                 localhomtable[i][j].nokori = 0;
 578                                 localhomtable[i][j].extended = -1;
 579                                 localhomtable[i][j].last = localhomtable[i]+j;
 580                                 localhomtable[i][j].korh = 'h';
 581                         }
 582                         if( !specifictarget ) ilim--;
 583                 }
 584                 fprintf( stderr, "Loading 'hat3' ... " );
 585                 fflush( stderr );
 586                 prep = fopen( "hat3", "r" );
 587                 if( prep == NULL ) ErrorExit( "Make hat3." );
 588                 if( specifictarget ) readlocalhomtable2_target( prep, njob, localhomtable, kozoarivec, targetmap );
 589                 else readlocalhomtable2_half( prep, njob, localhomtable, kozoarivec );
 590                 fclose( prep );
 591 //              for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 592 //                      fprintf( stdout, "%d %d %d %d %d %d %d\n", i, j, localhomtable[i][j].opt, localhomtable[i][j].start1, localhomtable[i][j].end1, localhomtable[i][j].start2, localhomtable[i][j].end2 );
 593                 fprintf( stderr, "done.\n" );
 594                 fflush( stderr );
 595 #if 0
 596                 fprintf( stderr, "Extending localhom ... " );
 597                 extendlocalhom( njob, localhomtable );
 598                 fprintf( stderr, "done.\n" );
 599 #endif
 600                 nkozo = 0;
 601                 for( i=0; i<njob; i++ ) if( kozoarivec[i] ) nkozo++;
 602         }
 603
 604
 605 //              if( nlenmax > 30000 )
 606                 if( nlenmax > 50000 ) // version >= 6.823
 607                 {
 608 #if 0
 609                         if( constraint )
 610                         {
 611                                 fprintf( stderr, "\nnlenmax=%d, nagasugi!\n", nlenmax );
 612                                 exit( 1 );
 613                         }
 614                         if( nevermemsave )
 615                         {
 616                                 fprintf( stderr, "\nnevermemsave=1, nlenmax=%d, nagasugi!\n", nlenmax );
 617                                 exit( 1 );
 618                         }
 619 #endif
 620                         if( !constraint && !nevermemsave && alg != 'M' )
 621                         {
 622                                 fprintf( stderr, "\nnlenmax=%d, Switching to the memsave mode\n", nlenmax );
 623                                 alg = 'M';
 624                         }
 625                 }
 626
 627
 628         if( specificityconsideration ) calcmaxdistclass();
 629
 630         for( i=0; i<njob; i++ )
 631         {
 632                 res_g[i][0] = 0;
 633         }
 634
 635         identity = 1;
 636         for( i=0; i<njob; i++ )
 637         {
 638                 identity *= ( nlen[i] == nlen[0] );
 639         }
 640         if( !identity )
 641         {
 642                 fprintf( stderr, "Input pre-aligned data\n" );
 643                 exit( 1 );
 644         }
 645         constants( njob, seq_g );
 646
 647 #if 0
 648         fprintf( stderr, "penalty = %d\n", penalty );
 649         fprintf( stderr, "penalty_ex = %d\n", penalty_ex );
 650         fprintf( stderr, "offset = %d\n", offset );
 651 #endif
 652
 653         initSignalSM();
 654
 655         initFiles();
 656
 657 #if 0
 658         if( njob == 2 )
 659         {
 660                 writePre( njob, name, nlen, seq_g, 1 );
 661                 SHOWVERSION;
 662                 return( 0 );
 663         }
 664 #else
 665         if( njob == 2 )
 666         {
 667                 weight = 0;
 668                 niter = 1;
 669         }
 670 #endif
 671
 672         c = seqcheck( seq_g );
 673         if( c )
 674         {
 675                 fprintf( stderr, "Illegal character %c\n", c );
 676                 exit( 1 );
 677         }
 678         commongappick( njob, seq_g );
 679
 680         if( rnakozo && rnaprediction == 'm' )
 681         {
 682                 singlerna = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 683                 prep = fopen( "hat4", "r" );
 684                 if( prep == NULL ) ErrorExit( "Make hat4 using mccaskill." );
 685                 fprintf( stderr, "Loading 'hat4' ... " );
 686                 fflush( stderr );
 687                 for( i=0; i<njob; i++ )
 688                 {
 689                         gappick0( nogap1seq[0], seq_g[i] );
 690                         nogaplen = strlen( nogap1seq[0] );
 691                         singlerna[i] = (RNApair **)calloc( nogaplen+1, sizeof( RNApair * ) );
 692                         for( j=0; j<nogaplen; j++ )
 693                         {
 694                                 singlerna[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
 695                                 singlerna[i][j][0].bestpos = -1;
 696                                 singlerna[i][j][0].bestscore = -1.0;
 697                         }
 698                         singlerna[i][nogaplen] = NULL;
 699                         readmccaskill( prep, singlerna[i], nogaplen );
 700                 }
 701                 fclose( prep );
 702                 fprintf( stderr, "\ndone.\n" );
 703                 fflush( stderr );
 704         }
 705         else
 706                 singlerna = NULL;
 707
 708
 709
 710         if( utree )
 711         {
 712                 prep = fopen( "hat2", "r" );
 713                 if( !prep ) ErrorExit( "Make hat2." );
 714                 readhat2_pointer( prep, njob, name, eff );
 715                 fclose( prep );
 716 #if 0
 717                 fprintf( stderr, "eff = \n" );
 718                 for( i=0; i<njob-1; i++ )
 719                 {
 720                         for( j=i+1; j<njob; j++ )
 721                         {
 722                                 fprintf( stderr, "%d-%d,  %f\n", i, j, eff[i][j] );
 723                         }
 724                         fprintf( stderr, "\n" );
 725                 }
 726 #endif
 727                 if( intree )
 728                 {
 729                         veryfastsupg_double_loadtree( njob, eff, topol, len, name );
 730 #if 0
 731                         fprintf( stderr, "eff = \n" );
 732                         for( i=0; i<njob-1; i++ )
 733                         {
 734                                 for( j=i+1; j<njob; j++ )
 735                                 {
 736                                         fprintf( stderr, "%d-%d,  %f\n", i, j, eff[i][j] );
 737                                 }
 738                                 fprintf( stderr, "\n" );
 739                         }
 740 exit( 1 );
 741 #endif
 742                 }
 743                 else if( intop ) // v6.528 deha if( intop ) dattanode intree ga mukou datta.
 744                 {
 745                         fprintf( stderr, "--topin has been disabled\n" );
 746                         exit( 1 );
 747 //                      veryfastsupg_double_loadtop( njob, eff, topol, len );
 748                 }
 749                 else if( subalignment )
 750                         fixed_supg_double_treeout_constrained( njob, eff, topol, len, name, nsubalignments, subtable );
 751                 else if( treemethod == 'X' || treemethod == 'E' || treemethod == 'q' )
 752 //                      veryfastsupg_double_outtree( njob, eff, topol, len, name );
 753                         fixed_musclesupg_double_treeout( njob, eff, topol, len, name );
 754                 else if( treemethod == 'n' )
 755                         nj( njob, eff, topol, len );
 756                 else if( treemethod == 's' )
 757                         spg( njob, eff, topol, len );
 758                 else if( treemethod == 'p' )
 759                         upg2( njob, eff, topol, len );
 760                 else ErrorExit( "Incorrect treemethod.\n" );
 761         }
 762 #if DEBUG
 763         printf( "utree = %d\n", utree );
 764 #endif
 765
 766         if( autosubalignment > 0.0 && subalignment == 0 )
 767         {
 768 //              reporterr( "Computing skipthisbranch..\n" );
 769                 insubtable = AllocateIntVec( njob );
 770                 preservegaps = AllocateIntVec( njob );
 771                 subtable = calloc( 1, sizeof( char * ) );
 772                 subtable[0] = NULL; // for FreeIntMtx
 773                 for( i=0; i<njob; i++ ) insubtable[i] = 0;
 774                 for( i=0; i<njob; i++ ) preservegaps[i] = 0; // tsukawanaikamo
 775                 if( generatesubalignmentstable( njob, &subtable, &nsubalignments, &maxmem, topol, len, autosubalignment ) ) // subtable ha allocate sareru.
 776                 {
 777                         reporterr( "################################################################################################ \n" );
 778                         reporterr( "#\n" );
 779                         reporterr( "# WARNING: Iterative refinment was not done because you gave a large --fix value (%6.3f).\n", autosubalignment );
 780                         reporterr( "#\n" );
 781                         reporterr( "################################################################################################ \n" );
 782                         writePre( njob, name, nlen, seq_g, 1 );
 783
 784                         FreeCharMtx( seq_g_bk );
 785                         FreeIntCub( topol );
 786                         FreeDoubleMtx( len );
 787                         FreeDoubleMtx( eff );
 788                         FreeIntMtx( skipthisbranch );
 789                         FreeIntMtx( subtable );
 790                         free( preservegaps );
 791                         free( insubtable );
 792                         SHOWVERSION;
 793                         return( 0 );
 794                 }
 795 //              subtable = AllocateIntMtx( nsubalignments, maxmem+1 );
 796                 fprintf( stderr, "nsubalignments = %d, maxmem = %d\n", nsubalignments, maxmem );
 797                 subalnpt = AllocateCharCub( nsubalignments, maxmem, 0 );
 798 #if 0
 799                 for( i=0; i<nsubalignments; i++ )
 800                 {
 801                         reporterr( "subalignment %d\n", i );
 802                         for( j=0; subtable[i][j]!=-1; j++ )
 803                         {
 804                                 reporterr( "%5d", subtable[i][j] );
 805                         }
 806                         reporterr( "\n" );
 807                 }
 808 #endif
 809 #if 0 // wakaran
 810                 for( i=0; i<nsubalignments; i++ ) for( j=0; j<insubtable[i]; j++ )
 811                 {
 812                         if( subtable[i][j] < 0 )
 813                         {
 814                                 fprintf( stderr, "Not supported in the iterative refinmenment mode.\n" );
 815                                 fprintf( stderr, "Please use a positive number to represent a sequence.\n" );
 816                         }
 817                 }
 818 #endif
 819 //              reporterr( "done.\n" );
 820         }
 821
 822
 823         orderfp = fopen( "order", "w" );
 824         if( !orderfp )
 825         {
 826                 fprintf( stderr, "Cannot open 'order'\n" );
 827                 exit( 1 );
 828         }
 829         for( i=0; (j=topol[njob-2][0][i])!=-1; i++ )
 830         {
 831                 fprintf( orderfp, "%d\n", j );
 832         }
 833         for( i=0; (j=topol[njob-2][1][i])!=-1; i++ )
 834         {
 835                 fprintf( orderfp, "%d\n", j );
 836         }
 837         fclose( orderfp );
 838
 839
 840
 841         fprintf( stderr, "\n" );
 842         if( ( !utree && kobetsubunkatsu ) || constraint || !bunkatsu )
 843         {
 844                 nseg = 0;
 845                 anchors[0] =0;
 846                 anchors[1] =strlen( seq_g[0] );
 847                 nseg += 2;
 848         }
 849         else
 850         {
 851                 nseg = searchAnchors( njob, seq_g, segment );
 852 #if 0
 853                 fprintf( stderr, "### nseg = %d\n", nseg );
 854                 fprintf( stderr, "### seq_g[0] = %s\n", seq_g[0] );
 855                 fprintf( stderr, "### nlenmax = %d\n", nlenmax );
 856                 fprintf( stderr, "### len = %d\n", strlen( seq_g[0] ) );
 857
 858 #endif
 859
 860                 anchors[0] = 0;
 861                 for( i=0; i<nseg; i++ ) anchors[i+1] = segment[i].center;
 862                 anchors[nseg+1] = strlen( seq_g[0] );
 863                 nseg +=2;
 864
 865 #if 0
 866                 for( i=0; i<nseg; i++ )
 867                         fprintf( stderr, "anchor[%d] = %d\n", i, anchors[i] );
 868 #endif
 869         }
 870
 871
 872 //--------------- kokokara ----
 873         if( subalignment || autosubalignment )
 874         {
 875                 for( i=0; i<nsubalignments; i++ )
 876                 {
 877                         fprintf( stderr, "\nChecking subalignment %d:\n", i+1 );
 878                         alignmentlength = strlen( seq[subtable[i][0]] );
 879                         for( j=0; subtable[i][j]!=-1; j++ )
 880                                 fprintf( stderr, " %d ", subtable[i][j]+1 );
 881 //                              fprintf( stderr, " ##### %d-%d. %-30.30s\n", i, subtable[i][j]+1, name[subtable[i][j]]+1 );
 882                         fprintf( stderr, "\n" );
 883                         for( j=0; subtable[i][j]!=-1; j++ )
 884                         {
 885                                 if( subtable[i][j] >= njob )
 886                                 {
 887                                         fprintf( stderr, "No such sequence, %d.\n", subtable[i][j]+1 );
 888                                         exit( 1 );
 889                                 }
 890                                 if( alignmentlength != strlen( seq[subtable[i][j]] ) )
 891                                 {
 892                                         fprintf( stderr, "\n" );
 893                                         fprintf( stderr, "###############################################################################\n" );
 894                                         fprintf( stderr, "# ERROR!\n" );
 895                                         fprintf( stderr, "# Subalignment %d must be aligned.\n", i+1 );
 896                                         fprintf( stderr, "# Please check the alignment lengths of following sequences.\n" );
 897                                         fprintf( stderr, "#\n" );
 898                                         fprintf( stderr, "# %d. %-10.10s -> %d letters (including gaps)\n", subtable[i][0]+1, name[subtable[i][0]]+1, alignmentlength );
 899                                         fprintf( stderr, "# %d. %-10.10s -> %d letters (including gaps)\n", subtable[i][j]+1, name[subtable[i][j]]+1, (int)strlen( seq[subtable[i][j]] ) );
 900                                         fprintf( stderr, "#\n" );
 901                                         fprintf( stderr, "# See http://mafft.cbrc.jp/alignment/software/merge.html for details.\n" );
 902                                         if( subalignmentoffset )
 903                                         {
 904                                                 fprintf( stderr, "#\n" );
 905                                                 fprintf( stderr, "# You specified seed alignment(s) consisting of %d sequences.\n", subalignmentoffset );
 906                                                 fprintf( stderr, "# In this case, the rule of numbering is:\n" );
 907                                                 fprintf( stderr, "#   The aligned seed sequences are numbered as 1 .. %d\n", subalignmentoffset );
 908                                                 fprintf( stderr, "#   The input sequences to be aligned are numbered as %d .. %d\n", subalignmentoffset+1, subalignmentoffset+njob );
 909                                         }
 910                                         fprintf( stderr, "###############################################################################\n" );
 911                                         fprintf( stderr, "\n" );
 912                                         exit( 1 );
 913                                 }
 914                                 insubtable[subtable[i][j]] = 1;
 915                         }
 916                         for( j=0; j<njob-1; j++ )
 917                         {
 918 #if 0
 919                                 int k;
 920                                 reporterr( "####  STEP%d\n", j );
 921                                 for( k=0; topol[j][0][k]!=-1; k++ ) reporterr( "%3d ", topol[j][0][k] );
 922                                 reporterr( "len=%f\n", len[j][0] );
 923                                 for( k=0; topol[j][1][k]!=-1; k++ ) reporterr( "%3d ", topol[j][1][k] );
 924                                 reporterr( "len=%f\n", len[j][1] );
 925                                 reporterr( "\n" );
 926 #endif
 927                                 if( includemember( topol[j][0], subtable[i] ) && !samemember( topol[j][0], subtable[i] ) )
 928                                 {
 929                                         skipthisbranch[j][0] = 1;
 930 //                                      reporterr( "SKIP 0 !!!!!!\n" );
 931                                 }
 932                                 if( includemember( topol[j][1], subtable[i] ) && !samemember( topol[j][1], subtable[i] ) )
 933                                 {
 934                                         skipthisbranch[j][1] = 1;
 935 //                                      reporterr( "SKIP 1 !!!!!!\n" );
 936                                 }
 937                         }
 938                         foundthebranch = 0;
 939                         for( j=0; j<njob-1; j++ )
 940                         {
 941                                 if( samemember( topol[j][0], subtable[i] ) || samemember( topol[j][1], subtable[i] ) )
 942                                 {
 943                                         foundthebranch = 1;
 944                                         fprintf( stderr, " -> OK\n" );
 945                                         break;
 946                                 }
 947                         }
 948                         if( !foundthebranch )
 949                         {
 950                                 system( "cp infile.tree GuideTree" ); // tekitou
 951                                 fprintf( stderr, "\n" );
 952                                 fprintf( stderr, "###############################################################################\n" );
 953                                 fprintf( stderr, "# ERROR!\n" );
 954                                 fprintf( stderr, "# Subalignment %d does not seem to form a monophyletic cluster\n", i+1 );
 955                                 fprintf( stderr, "# in the guide tree ('GuideTree' in this directory) internally computed.\n" );
 956                                 fprintf( stderr, "# If you really want to use this subalignment, pelase give a tree with --treein \n" );
 957                                 fprintf( stderr, "# http://mafft.cbrc.jp/alignment/software/treein.html\n" );
 958                                 fprintf( stderr, "# http://mafft.cbrc.jp/alignment/software/merge.html\n" );
 959                                 if( subalignmentoffset )
 960                                 {
 961                                         fprintf( stderr, "#\n" );
 962                                         fprintf( stderr, "# You specified seed alignment(s) consisting of %d sequences.\n", subalignmentoffset );
 963                                         fprintf( stderr, "# In this case, the rule of numbering is:\n" );
 964                                         fprintf( stderr, "#   The aligned seed sequences are numbered as 1 .. %d\n", subalignmentoffset );
 965                                         fprintf( stderr, "#   The input sequences to be aligned are numbered as %d .. %d\n", subalignmentoffset+1, subalignmentoffset+njob );
 966                                 }
 967                                 fprintf( stderr, "############################################################################### \n" );
 968                                 fprintf( stderr, "\n" );
 969                                 exit( 1 );
 970                         }
 971 //                      commongappick( seq[subtable[i]], subalignment[i] ); // irukamo
 972                 }
 973 #if 0
 974                 for( i=0; i<njob-1; i++ )
 975                 {
 976                         fprintf( stderr, "STEP %d\n", i+1 );
 977                         fprintf( stderr, "group1 = " );
 978                         for( j=0; topol[i][0][j] != -1; j++ )
 979                                 fprintf( stderr, "%d ", topol[i][0][j]+1 );
 980                         fprintf( stderr, "\n" );
 981                         fprintf( stderr, "SKIP -> %d\n\n", skipthisbranch[i][0] );
 982                         fprintf( stderr, "group2 = " );
 983                         for( j=0; topol[i][1][j] != -1; j++ )
 984                                 fprintf( stderr, "%d ", topol[i][1][j]+1 );
 985                         fprintf( stderr, "\n" );
 986                         fprintf( stderr, "SKIP -> %d\n\n", skipthisbranch[i][1] );
 987                 }
 988 #endif
 989
 990                 for( i=0; i<njob; i++ )
 991                 {
 992                         if( insubtable[i] ) strcpy( bseq[i], seq[i] );
 993                         else gappick0( bseq[i], seq[i] );
 994                 }
 995
 996                 for( i=0; i<nsubalignments; i++ )
 997                 {
 998                         for( j=0; subtable[i][j]!=-1; j++ ) subalnpt[i][j] = bseq[subtable[i][j]];
 999                         commongappick( j, subalnpt[i] );
1000                 }
1001
1002                 FreeIntMtx( subtable );
1003                 free( insubtable );
1004                 for( i=0; i<nsubalignments; i++ ) free( subalnpt[i] );
1005                 free( subalnpt );
1006                 free( preservegaps );
1007         }
1008 //--------------- kokomade ----
1009
1010
1011
1012
1013         for( i=0; i<njob; i++ ) res_g[i][0] = 0;
1014
1015         for( iseg=0; iseg<nseg-1; iseg++ )
1016         {
1017                 int tmplen = anchors[iseg+1]-anchors[iseg];
1018                 int pos = strlen( res_g[0] );
1019                 for( j=0; j<njob; j++ )
1020                 {
1021                         strncpy( seq[j], seq_g[j], tmplen );
1022                         seq[j][tmplen]= 0;
1023                         seq_g[j] += tmplen;
1024
1025                 }
1026                 fprintf( stderr, "Segment %3d/%3d %4d-%4d\n", iseg+1, nseg-1, pos+1, pos+1+tmplen );
1027                 fflush( stderr );
1028                 fprintf( trap_g, "Segment %3d/%3d %4d-%4d\n", iseg+1, nseg-1, pos+1, pos+1+tmplen );
1029
1030                 cut = ocut;
1031                 returnvalue = TreeDependentIteration( njob, name, nlen, seq, bseq, topol, len, eff, skipthisbranch, alloclen, localhomtable, singlerna, nkozo, kozoarivec, ntarget, targetmap, targetmapr );
1032
1033                 for( i=0; i<njob; i++ )
1034                         strcat( res_g[i], bseq[i] );
1035         }
1036         FreeCharMtx( seq_g_bk );
1037         FreeIntCub( topol );
1038         FreeDoubleMtx( len );
1039         FreeDoubleMtx( eff );
1040         FreeIntMtx( skipthisbranch );
1041         free( kozoarivec );
1042         if( constraint )
1043         {
1044                 if( specifictarget ) FreeLocalHomTable_part( localhomtable, ntarget, njob );
1045                 else FreeLocalHomTable_half( localhomtable, njob );
1046         }
1047         free( targetmap );
1048         free( targetmapr );
1049         if( rnakozo && rnaprediction == 'm' )
1050         {
1051                 if( singlerna ) // nen no tame
1052                 {
1053                         for( i=0; i<njob; i++ )
1054                         {
1055                                 for( j=0; singlerna[i][j]!=NULL; j++ )
1056                                 {
1057                                         if( singlerna[i][j] ) free( singlerna[i][j] );
1058                                 }
1059                                 if( singlerna[i] ) free( singlerna[i] );
1060                         }
1061                         free( singlerna );
1062                         singlerna = NULL;
1063                 }
1064         }
1065
1066 #if 0
1067         Write( stdout, njob, name, nlen, bseq );
1068 #endif
1069
1070         fprintf( stderr, "done\n" );
1071         fprintf( trap_g, "done\n" );
1072         fclose( trap_g );
1073         freeconstants();
1074
1075
1076         devide = 0;
1077         writePre( njob, name, nlen, res_g, 1 );
1078 #if 0
1079         writeData( stdout, njob, name, nlen, res_g, 1 );
1080 #endif
1081
1082
1083         if( spscoreout ) reporterr( "Unweighted sum-of-pairs score = %10.5f\n", sumofpairsscore( njob, res_g ) );
1084
1085         SHOWVERSION;
1086         return( 0 );
1087 }
1088
1089 #if 0
1090 signed int main( int argc, char *argv[] )
1091 {
1092         int i, nlen[M];
1093         char b[B];
1094         char a[] = "=";
1095         int value;
1096
1097         gets( b ); njob = atoi( b );
1098
1099 /*
1100         scoremtx = 0;
1101         if( strstr( b, "ayhoff" ) ) scoremtx = 1;
1102         else if( strstr( b, "dna" ) || strstr( b, "DNA" ) ) scoremtx = -1;
1103         else if( strstr( b, "M-Y" ) || strstr( b, "iyata" ) ) scoremtx = 2;
1104         else scoremtx = 0;
1105 */
1106         if( strstr( b, "constraint" ) ) cnst = 1;
1107
1108         nlenmax = 0;
1109         i = 0;
1110         while( i<njob )
1111         {
1112                 gets( b );
1113                 if( !strncmp( b, a, 1 ) )
1114                 {
1115                         gets( b ); nlen[i] = atoi( b );
1116                         if( nlen[i] > nlenmax ) nlenmax = nlen[i];
1117                         i++;
1118                 }
1119         }
1120         if( nlenmax > N || njob > M )
1121         {
1122                 fprintf( stderr, "ERROR in main\n" );
1123                 exit( 1 );
1124         }
1125         /*
1126         nlenmax = Na;
1127         */
1128         rewind( stdin );
1129         value = main1( nlen, argc, argv );
1130         exit( 0 );
1131 }
1132 #endif