binaries/src/mafft/core/getlag.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5
   6 void arguments( int argc, char *argv[] )
   7 {
   8     int c;
   9
  10         calledByXced = 0;
  11         devide = 0;
  12         use_fft = 0;
  13         fftscore = 1;
  14         fftRepeatStop = 0;
  15         fftNoAnchStop = 0;
  16     weight = 3;
  17     utree = 1;
  18         tbutree = 1;
  19     refine = 0;
  20     check = 1;
  21     cut = 0.0;
  22     disp = 0;
  23     outgap = 1;
  24     alg = 'C';
  25     mix = 0;
  26         tbitr = 0;
  27         scmtd = 5;
  28         tbweight = 0;
  29         tbrweight = 3;
  30         checkC = 0;
  31         treemethod = 'x';
  32         contin = 0;
  33         ppenalty = NOTSPECIFIED;
  34         ppenalty_ex = NOTSPECIFIED;
  35         poffset = NOTSPECIFIED;
  36         kimuraR = NOTSPECIFIED;
  37         pamN = NOTSPECIFIED;
  38         geta2 = GETA2;
  39         scoremtx = NOTSPECIFIED;
  40
  41     while( --argc > 0 && (*++argv)[0] == '-' )
  42         {
  43         while ( (c = *++argv[0]) )
  44                 {
  45             switch( c )
  46             {
  47                                 case 'f':
  48                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
  49                                         fprintf( stderr, "ppenalty = %d\n", ppenalty );
  50                                         --argc;
  51                                         goto nextoption;
  52                                 case 'g':
  53                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
  54                                         fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
  55                                         --argc;
  56                                         goto nextoption;
  57                                 case 'h':
  58                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
  59                                         fprintf( stderr, "poffset = %d\n", poffset );
  60                                         --argc;
  61                                         goto nextoption;
  62                                 case 'D':
  63                                         scoremtx = -1;
  64                                         break;
  65                                 case 'P':
  66                                         scoremtx = 0;
  67                                         break;
  68                                 case 'i':
  69                                         contin = 1;
  70                                         break;
  71                                 case 'e':
  72                                         fftscore = 0;
  73                                         break;
  74                                 case 'O':
  75                                         fftNoAnchStop = 1;
  76                                         break;
  77                                 case 'R':
  78                                         fftRepeatStop = 1;
  79                                         break;
  80                                 case 'Q':
  81                                         calledByXced = 1;
  82                                         break;
  83                                 case 's':
  84                                         treemethod = 's';
  85                                         break;
  86                                 case 'x':
  87                                         treemethod = 'x';
  88                                         break;
  89                                 case 'p':
  90                                         treemethod = 'p';
  91                                         break;
  92                                 case 'a':
  93                                         alg = 'a';
  94                                         break;
  95                                 case 'A':
  96                                         alg = 'A';
  97                                         break;
  98                                 case 'S':
  99                                         alg = 'S';
 100                                         break;
 101                                 case 'C':
 102                                         alg = 'C';
 103                                         break;
 104                                 case 'F':
 105                                         use_fft = 1;
 106                                         break;
 107                                 case 'v':
 108                                         tbrweight = 3;
 109                                         break;
 110                                 case 'd':
 111                                         disp = 1;
 112                                         break;
 113                                 case 'o':
 114                                         outgap = 0;
 115                                         break;
 116 /* Modified 01/08/27, default: user tree */
 117                                 case 'J':
 118                                         tbutree = 0;
 119                                         break;
 120 /* modification end. */
 121                                 case 'Z':
 122                                         checkC = 1;
 123                                         break;
 124                 default:
 125                     fprintf( stderr, "illegal option %c\n", c );
 126                     argc = 0;
 127                     break;
 128             }
 129                 }
 130                 nextoption:
 131                         ;
 132         }
 133     if( argc == 1 )
 134     {
 135         cut = atof( (*argv) );
 136         argc--;
 137     }
 138     if( argc != 0 )
 139     {
 140         fprintf( stderr, "options: Check source file !\n" );
 141         exit( 1 );
 142     }
 143         if( tbitr == 1 && outgap == 0 )
 144         {
 145                 fprintf( stderr, "conflicting options : o, m or u\n" );
 146                 exit( 1 );
 147         }
 148         if( alg == 'C' && outgap == 0 )
 149         {
 150                 fprintf( stderr, "conflicting options : C, o\n" );
 151                 exit( 1 );
 152         }
 153         readOtherOptions( &ppid, &fftThreshold, &fftWinSize );
 154 }
 155
 156
 157 void treebase( char name[M][B], int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, double **mtx, int ***topol, double **len, double **eff, int alloclen )
 158 {
 159         int i, j, l;
 160         int clus1, clus2;
 161         int s1, s2, r1, r2;
 162         float pscore;
 163         static char *indication1, *indication2;
 164         static char name1[M][B], name2[M][B];
 165         static double **partialmtx = NULL;
 166         static int ***partialtopol = NULL;
 167         static double **partiallen = NULL;
 168         static double **partialeff = NULL;
 169         static double *effarr = NULL;
 170         static double *effarr1 = NULL;
 171         static double *effarr2 = NULL;
 172 #if 0
 173         char pair[njob][njob];
 174 #else
 175         static char **pair;
 176 #endif
 177         if( partialtopol == NULL )
 178         {
 179                 partialmtx = AllocateDoubleMtx( njob, njob );
 180                 partialtopol = AllocateIntCub( njob, 2, njob );
 181                 partialeff = AllocateDoubleMtx( njob, njob );
 182                 partiallen = AllocateDoubleMtx( njob, 2 );
 183                 effarr = AllocateDoubleVec( njob );
 184                 effarr1 = AllocateDoubleVec( njob );
 185                 effarr2 = AllocateDoubleVec( njob );
 186                 indication1 = AllocateCharVec( njob*3+100 );
 187                 indication2 = AllocateCharVec( njob*3+100 );
 188 #if 0
 189 #else
 190                 pair = AllocateCharMtx( njob, njob );
 191 #endif
 192         }
 193
 194         if( checkC )
 195                 for( i=0; i<njob; i++ ) fprintf( stderr, "eff in tb-%d %f\n", i, eff[i][i] );
 196
 197         for( i=0; i<njob; i++ ) effarr[i] = eff[i][i];
 198         for( i=0; i<njob; i++ ) strcpy( aseq[i], seq[i] );
 199
 200         if( checkC )
 201                 for( i=0; i<njob; i++ ) fprintf( stderr, "effarr for aseq-%d %f\n", i, effarr[i] );
 202
 203         writePre( njob, name, nlen, aseq, 0 );
 204
 205         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
 206         for( i=0; i<njob; i++ ) pair[i][i] = 1;
 207         for( l=0; l<njob-1; l++ )
 208         {
 209                 s1 = topol[l][0][0];
 210                 for( i=0; (r1=topol[l][0][i])>-1; i++ )
 211                         if( pair[s1][r1] != 1 ) exit( 1 );
 212                 s2 = topol[l][1][0];
 213                 for( i=0; (r2=topol[l][1][i])>-1; i++ )
 214                         if( pair[s2][r2] != 1 ) exit( 1 );
 215
 216                 clus1 = conjuction( pair, s1, aseq, mseq1, effarr1, effarr, name, name1, indication1 );
 217                 clus2 = conjuction( pair, s2, aseq, mseq2, effarr2, effarr, name, name2, indication2 );
 218                 fprintf( trap_g, "\nSTEP-%d\n", l );
 219                 fprintf( trap_g, "group1 = %s\n", indication1 );
 220                 fprintf( trap_g, "group2 = %s\n", indication2 );
 221
 222                 fprintf( stderr, "STEP %d /%d\n", l+1, njob-1 );
 223                 fprintf( stderr, "group1 = %.66s", indication1 );
 224                 if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
 225                 fprintf( stderr, "\n" );
 226                 fprintf( stderr, "group2 = %.66s", indication2 );
 227                 if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
 228                 fprintf( stderr, "\n" );
 229
 230                 if( checkC )
 231                         for( i=0; i<clus1; i++ ) fprintf( stderr, "STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
 232 /*
 233                 fprintf( stderr, "before align all\n" );
 234                 display( aseq, njob );
 235                 fprintf( stderr, "\n" );
 236                 fprintf( stderr, "before align 1 %s \n", indication1 );
 237                 display( mseq1, clus1 );
 238                 fprintf( stderr, "\n" );
 239                 fprintf( stderr, "before align 2 %s \n", indication2 );
 240                 display( mseq2, clus2 );
 241                 fprintf( stderr, "\n" );
 242 */
 243
 244                 pscore = Fgetlag( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 245
 246                 for( i=0; (r2=topol[l][1][i])>-1; i++ )
 247                 {
 248                         pair[s1][r2] = 1;
 249                         pair[s2][r2] = 0;
 250                 }
 251
 252                 writePre( njob, name, nlen, aseq, 0 );
 253
 254                 if( disp ) display( aseq, njob );
 255                 fprintf( stderr, "\n" );
 256
 257         }
 258 }
 259
 260 static void WriteOptions( FILE *fp )
 261 {
 262         fprintf( fp, "tree-base method\n" );
 263         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
 264         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
 265         if( tbitr || tbweight )
 266         {
 267                 fprintf( fp, "iterate at each step\n" );
 268                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
 269                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
 270                 if( tbweight ) fprintf( fp, "  weighted\n" );
 271                 fprintf( fp, "\n" );
 272         }
 273         if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
 274         else if( scoremtx ==  1 ) fprintf( fp, "Dayhoff( machigai ga aru )\n" );
 275         else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
 276         else if( scoremtx == -1 ) fprintf( fp, "DNA\n" );
 277
 278     if( scoremtx == 0 || scoremtx == -1 )
 279         fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 280     else
 281         fprintf( fp, "Gap Penalty = %+5.2f\n", (double)ppenalty/1000 );
 282
 283         if( alg == 'a' )
 284                 fprintf( fp, "Algorithm A\n" );
 285         else if( alg == 'A' )
 286                 fprintf( fp, "Apgorithm A+\n" );
 287         else if( alg == 'S' )
 288                 fprintf( fp, "Apgorithm S\n" );
 289         else if( alg == 'C' )
 290                 fprintf( fp, "Apgorithm A+/C\n" );
 291         else
 292                 fprintf( fp, "Unknown algorithm\n" );
 293
 294         if( treemethod == 'x' )
 295                 fprintf( fp, "Tree = UPGMA (3).\n" );
 296         else if( treemethod == 's' )
 297                 fprintf( fp, "Tree = UPGMA (2).\n" );
 298         else if( treemethod == 'p' )
 299                 fprintf( fp, "Tree = UPGMA (1).\n" );
 300         else
 301                 fprintf( fp, "Unknown tree.\n" );
 302
 303     if( use_fft )
 304     {
 305         fprintf( fp, "FFT on\n" );
 306         if( scoremtx == -1 )
 307             fprintf( fp, "Basis : 4 nucleotides\n" );
 308         else
 309         {
 310             if( fftscore )
 311                 fprintf( fp, "Basis : Polarity and Volume\n" );
 312             else
 313                 fprintf( fp, "Basis : 20 amino acids\n" );
 314         }
 315         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
 316         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
 317     }
 318         else
 319         fprintf( fp, "FFT off\n" );
 320         fflush( fp );
 321 }
 322
 323
 324 int main( int argc, char *argv[] )
 325 {
 326         static int  nlen[M];
 327         static char name[M][B], **seq;
 328         static char **mseq1, **mseq2;
 329         static char **aseq;
 330         static char **bseq;
 331         static double **pscore;
 332         static double **eff;
 333         static double **node0, **node1;
 334         int i, j;
 335         static int ***topol;
 336         static double **len;
 337         FILE *prep;
 338         char c;
 339         int alloclen;
 340
 341         arguments( argc, argv );
 342         getnumlen( stdin );
 343         rewind( stdin );
 344
 345         seq = AllocateCharMtx( njob, nlenmax*5+1 );
 346         aseq = AllocateCharMtx( njob, nlenmax*5+1 );
 347         bseq = AllocateCharMtx( njob, nlenmax*5+1 );
 348         mseq1 = AllocateCharMtx( njob, 0 );
 349         mseq2 = AllocateCharMtx( njob, 0 );
 350         alloclen = nlenmax*5;
 351
 352         topol = AllocateIntCub( njob, 2, njob );
 353         len = AllocateDoubleMtx( njob, 2 );
 354         pscore = AllocateDoubleMtx( njob, njob );
 355         eff = AllocateDoubleMtx( njob, njob );
 356         node0 = AllocateDoubleMtx( njob, njob );
 357         node1 = AllocateDoubleMtx( njob, njob );
 358
 359 #if 0
 360         Read( name, nlen, seq );
 361 #else
 362         readData( stdin, name, nlen, seq );
 363 #endif
 364
 365         constants( njob, seq );
 366
 367 #if 0
 368         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
 369 #endif
 370
 371         initSignalSM();
 372
 373         initFiles();
 374
 375         WriteOptions( trap_g );
 376
 377         c = seqcheck( seq );
 378         if( c )
 379         {
 380                 fprintf( stderr, "Illeagal character %c\n", c );
 381                 exit( 1 );
 382         }
 383
 384         writePre( njob, name, nlen, seq, 0 );
 385
 386         if( tbutree == 0 )
 387         {
 388                 for( i=1; i<njob; i++ )
 389                 {
 390                         if( nlen[i] != nlen[0] )
 391                         {
 392                                 fprintf( stderr, "Input pre-aligned seqences or make hat2.\n" );
 393                                 exit( 1 );
 394                         }
 395                 }
 396                 for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 397                 {
 398                 /*
 399                         pscore[i][j] = (double)score_calc1( seq[i], seq[j] );
 400                 */
 401                         pscore[i][j] = (double)substitution_hosei( seq[i], seq[j] );
 402                 }
 403         }
 404         else
 405         {
 406                 fprintf( stderr, "Loading 'hat2' ... " );
 407                 prep = fopen( "hat2", "r" );
 408                 if( prep == NULL ) ErrorExit( "Make hat2." );
 409                 readhat2( prep, njob, name, pscore );
 410                 fclose( prep );
 411                 fprintf( stderr, "done.\n" );
 412
 413 #if 0
 414                 prep = fopen( "hat2_check", "w" );
 415                 WriteHat2( prep, njob, name, pscore );
 416                 fclose( prep );
 417 #endif
 418
 419         }
 420
 421         fprintf( stderr, "Constructing dendrogram ... " );
 422         if( treemethod == 'x' )
 423                 supg( njob, pscore, topol, len );
 424         else if( treemethod == 's' )
 425                 spg( njob, pscore, topol, len );
 426         else if( treemethod == 'p' )
 427                 upg2( njob, pscore, topol, len );
 428         else
 429                 ErrorExit( "Incorrect tree\n" );
 430         fprintf( stderr, "done.\n" );
 431
 432         countnode( njob, topol, node0 );
 433         if( tbrweight )
 434         {
 435                 weight = 3;
 436                 utree = 0; counteff( njob, topol, len, eff ); utree = 1;
 437         }
 438         else
 439         {
 440                 for( i=0; i<njob; i++ ) eff[i][i] = 1.0;
 441         }
 442
 443
 444         for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
 445
 446         treebase( name, nlen, bseq, aseq, mseq1, mseq2, pscore, topol, len, eff, alloclen );
 447
 448         fprintf( trap_g, "done\n" );
 449         fclose( trap_g );
 450
 451         writePre( njob, name, nlen, aseq, !contin );
 452         writeData( stdout, njob, name, nlen, aseq );
 453 #if IODEBUG
 454         fprintf( stderr, "OSHIMAI\n" );
 455 #endif
 456         SHOWVERSION;
 457         return( 0 );
 458 }