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, 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, **name2;
 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                 name1 = AllocateCharMtx( njob, B+1 );
 189                 name2 = AllocateCharMtx( njob, B+1 );
 190 #if 0
 191 #else
 192                 pair = AllocateCharMtx( njob, njob );
 193 #endif
 194         }
 195
 196         if( checkC )
 197                 for( i=0; i<njob; i++ ) fprintf( stderr, "eff in tb-%d %f\n", i, eff[i][i] );
 198
 199         for( i=0; i<njob; i++ ) effarr[i] = eff[i][i];
 200         for( i=0; i<njob; i++ ) strcpy( aseq[i], seq[i] );
 201
 202         if( checkC )
 203                 for( i=0; i<njob; i++ ) fprintf( stderr, "effarr for aseq-%d %f\n", i, effarr[i] );
 204
 205         writePre( njob, name, nlen, aseq, 0 );
 206
 207         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
 208         for( i=0; i<njob; i++ ) pair[i][i] = 1;
 209         for( l=0; l<njob-1; l++ )
 210         {
 211                 s1 = topol[l][0][0];
 212                 for( i=0; (r1=topol[l][0][i])>-1; i++ )
 213                         if( pair[s1][r1] != 1 ) exit( 1 );
 214                 s2 = topol[l][1][0];
 215                 for( i=0; (r2=topol[l][1][i])>-1; i++ )
 216                         if( pair[s2][r2] != 1 ) exit( 1 );
 217
 218                 clus1 = conjuction( pair, s1, aseq, mseq1, effarr1, effarr, name, name1, indication1 );
 219                 clus2 = conjuction( pair, s2, aseq, mseq2, effarr2, effarr, name, name2, indication2 );
 220                 fprintf( trap_g, "\nSTEP-%d\n", l );
 221                 fprintf( trap_g, "group1 = %s\n", indication1 );
 222                 fprintf( trap_g, "group2 = %s\n", indication2 );
 223
 224                 fprintf( stderr, "STEP %d /%d\n", l+1, njob-1 );
 225                 fprintf( stderr, "group1 = %.66s", indication1 );
 226                 if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
 227                 fprintf( stderr, "\n" );
 228                 fprintf( stderr, "group2 = %.66s", indication2 );
 229                 if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
 230                 fprintf( stderr, "\n" );
 231
 232                 if( checkC )
 233                         for( i=0; i<clus1; i++ ) fprintf( stderr, "STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
 234 /*
 235                 fprintf( stderr, "before align all\n" );
 236                 display( aseq, njob );
 237                 fprintf( stderr, "\n" );
 238                 fprintf( stderr, "before align 1 %s \n", indication1 );
 239                 display( mseq1, clus1 );
 240                 fprintf( stderr, "\n" );
 241                 fprintf( stderr, "before align 2 %s \n", indication2 );
 242                 display( mseq2, clus2 );
 243                 fprintf( stderr, "\n" );
 244 */
 245
 246                 pscore = Fgetlag( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 247
 248                 for( i=0; (r2=topol[l][1][i])>-1; i++ )
 249                 {
 250                         pair[s1][r2] = 1;
 251                         pair[s2][r2] = 0;
 252                 }
 253
 254                 writePre( njob, name, nlen, aseq, 0 );
 255
 256                 if( disp ) display( aseq, njob );
 257                 fprintf( stderr, "\n" );
 258
 259         }
 260 }
 261
 262 static void WriteOptions( FILE *fp )
 263 {
 264         fprintf( fp, "tree-base method\n" );
 265         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
 266         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
 267         if( tbitr || tbweight )
 268         {
 269                 fprintf( fp, "iterate at each step\n" );
 270                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
 271                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
 272                 if( tbweight ) fprintf( fp, "  weighted\n" );
 273                 fprintf( fp, "\n" );
 274         }
 275         if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
 276         else if( scoremtx ==  1 ) fprintf( fp, "Dayhoff( machigai ga aru )\n" );
 277         else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
 278         else if( scoremtx == -1 ) fprintf( fp, "DNA\n" );
 279
 280     if( scoremtx == 0 || scoremtx == -1 )
 281         fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 282     else
 283         fprintf( fp, "Gap Penalty = %+5.2f\n", (double)ppenalty/1000 );
 284
 285         if( alg == 'a' )
 286                 fprintf( fp, "Algorithm A\n" );
 287         else if( alg == 'A' )
 288                 fprintf( fp, "Apgorithm A+\n" );
 289         else if( alg == 'S' )
 290                 fprintf( fp, "Apgorithm S\n" );
 291         else if( alg == 'C' )
 292                 fprintf( fp, "Apgorithm A+/C\n" );
 293         else
 294                 fprintf( fp, "Unknown algorithm\n" );
 295
 296         if( treemethod == 'x' )
 297                 fprintf( fp, "Tree = UPGMA (3).\n" );
 298         else if( treemethod == 's' )
 299                 fprintf( fp, "Tree = UPGMA (2).\n" );
 300         else if( treemethod == 'p' )
 301                 fprintf( fp, "Tree = UPGMA (1).\n" );
 302         else
 303                 fprintf( fp, "Unknown tree.\n" );
 304
 305     if( use_fft )
 306     {
 307         fprintf( fp, "FFT on\n" );
 308         if( scoremtx == -1 )
 309             fprintf( fp, "Basis : 4 nucleotides\n" );
 310         else
 311         {
 312             if( fftscore )
 313                 fprintf( fp, "Basis : Polarity and Volume\n" );
 314             else
 315                 fprintf( fp, "Basis : 20 amino acids\n" );
 316         }
 317         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
 318         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
 319     }
 320         else
 321         fprintf( fp, "FFT off\n" );
 322         fflush( fp );
 323 }
 324
 325
 326 int main( int argc, char *argv[] )
 327 {
 328         static int  nlen[M];
 329         static char **name, **seq;
 330         static char **mseq1, **mseq2;
 331         static char **aseq;
 332         static char **bseq;
 333         static double **pscore;
 334         static double **eff;
 335         static double **node0, **node1;
 336         int i, j;
 337         static int ***topol;
 338         static double **len;
 339         FILE *prep;
 340         char c;
 341         int alloclen;
 342
 343         arguments( argc, argv );
 344         getnumlen( stdin );
 345         rewind( stdin );
 346
 347         name = AllocateCharMtx( njob, B+1 );
 348         seq = AllocateCharMtx( njob, nlenmax*5+1 );
 349         aseq = AllocateCharMtx( njob, nlenmax*5+1 );
 350         bseq = AllocateCharMtx( njob, nlenmax*5+1 );
 351         mseq1 = AllocateCharMtx( njob, 0 );
 352         mseq2 = AllocateCharMtx( njob, 0 );
 353         alloclen = nlenmax*5;
 354
 355         topol = AllocateIntCub( njob, 2, njob );
 356         len = AllocateDoubleMtx( njob, 2 );
 357         pscore = AllocateDoubleMtx( njob, njob );
 358         eff = AllocateDoubleMtx( njob, njob );
 359         node0 = AllocateDoubleMtx( njob, njob );
 360         node1 = AllocateDoubleMtx( njob, njob );
 361
 362 #if 0
 363         Read( name, nlen, seq );
 364 #else
 365         readData_pointer( stdin, name, nlen, seq );
 366 #endif
 367
 368         constants( njob, seq );
 369
 370 #if 0
 371         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
 372 #endif
 373
 374         initSignalSM();
 375
 376         initFiles();
 377
 378         WriteOptions( trap_g );
 379
 380         c = seqcheck( seq );
 381         if( c )
 382         {
 383                 fprintf( stderr, "Illeagal character %c\n", c );
 384                 exit( 1 );
 385         }
 386
 387         writePre( njob, name, nlen, seq, 0 );
 388
 389         if( tbutree == 0 )
 390         {
 391                 for( i=1; i<njob; i++ )
 392                 {
 393                         if( nlen[i] != nlen[0] )
 394                         {
 395                                 fprintf( stderr, "Input pre-aligned seqences or make hat2.\n" );
 396                                 exit( 1 );
 397                         }
 398                 }
 399                 for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 400                 {
 401                 /*
 402                         pscore[i][j] = (double)score_calc1( seq[i], seq[j] );
 403                 */
 404                         pscore[i][j] = (double)substitution_hosei( seq[i], seq[j] );
 405                 }
 406         }
 407         else
 408         {
 409                 fprintf( stderr, "Loading 'hat2' ... " );
 410                 prep = fopen( "hat2", "r" );
 411                 if( prep == NULL ) ErrorExit( "Make hat2." );
 412                 readhat2_pointer( prep, njob, name, pscore );
 413                 fclose( prep );
 414                 fprintf( stderr, "done.\n" );
 415
 416 #if 0
 417                 prep = fopen( "hat2_check", "w" );
 418                 WriteHat2( prep, njob, name, pscore );
 419                 fclose( prep );
 420 #endif
 421
 422         }
 423
 424         fprintf( stderr, "Constructing dendrogram ... " );
 425         if( treemethod == 'x' )
 426                 supg( njob, pscore, topol, len );
 427         else if( treemethod == 's' )
 428                 spg( njob, pscore, topol, len );
 429         else if( treemethod == 'p' )
 430                 upg2( njob, pscore, topol, len );
 431         else
 432                 ErrorExit( "Incorrect tree\n" );
 433         fprintf( stderr, "done.\n" );
 434
 435         countnode( njob, topol, node0 );
 436         if( tbrweight )
 437         {
 438                 weight = 3;
 439                 utree = 0; counteff( njob, topol, len, eff ); utree = 1;
 440         }
 441         else
 442         {
 443                 for( i=0; i<njob; i++ ) eff[i][i] = 1.0;
 444         }
 445
 446
 447         for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
 448
 449         treebase( name, nlen, bseq, aseq, mseq1, mseq2, pscore, topol, len, eff, alloclen );
 450
 451         fprintf( trap_g, "done\n" );
 452         fclose( trap_g );
 453
 454         writePre( njob, name, nlen, aseq, !contin );
 455         writeData_pointer( stdout, njob, name, nlen, aseq );
 456 #if IODEBUG
 457         fprintf( stderr, "OSHIMAI\n" );
 458 #endif
 459         SHOWVERSION;
 460         return( 0 );
 461 }