binaries/src/mafft/core/mafft-profile.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4
   5 #if DEBUG
   6 #include <time.h>
   7 #include <sys/time.h>
   8 #include <sys/resource.h>
   9 double getrusage_sec()
  10 {
  11     struct rusage t;
  12     struct timeval tv;
  13     getrusage(RUSAGE_SELF, &t);
  14     tv = t.ru_utime;
  15     return tv.tv_sec + (double)tv.tv_usec*1e-6;
  16 }
  17 #endif
  18
  19
  20 int intcmp( int *str1, int *str2 )
  21 {
  22         while( *str1 != -1 && *str2 != -1 )
  23                 if( *str1++ != *str2++ ) return( 1 );
  24         if( *str1 != *str2 ) return( 1 );
  25         return( 0 );
  26 }
  27
  28 char **arguments( int argc, char *argv[] )
  29 {
  30     int c = 0;
  31
  32         fmodel = 0;
  33         nblosum = 62;
  34         calledByXced = 0;
  35         devide = 0;
  36         fftscore = 1;
  37         use_fft = 1;
  38         alg = 'A';
  39     weight = 0;
  40     utree = 1;
  41         tbutree = 0;
  42     refine = 0;
  43     check = 1;
  44     cut = 0.0;
  45     disp = 0;
  46     outgap = 0;
  47     mix = 0;
  48         tbitr = 0;
  49         scmtd = 5;
  50         tbweight = 0;
  51         tbrweight = 3;
  52         checkC = 0;
  53     scoremtx = 1;
  54         dorp = NOTSPECIFIED;
  55         ppenalty = NOTSPECIFIED;
  56         ppenalty_ex = NOTSPECIFIED;
  57         poffset = 0; // chokusetsu yobareru kara
  58         kimuraR = NOTSPECIFIED;
  59         pamN = NOTSPECIFIED;
  60         fftWinSize = NOTSPECIFIED;
  61         fftThreshold = NOTSPECIFIED;
  62         TMorJTT = JTT;
  63         treemethod = 'x';
  64
  65
  66     while( --argc > 0 && (*++argv)[0] == '-' )
  67         {
  68         while ( ( c = *++argv[0] ) )
  69                 {
  70             switch( c )
  71             {
  72                                 case 'P':
  73                                         dorp = 'p';
  74                                         break;
  75                                 case 'D':
  76                                         dorp = 'd';
  77                                         break;
  78                                 case 'F':
  79                                         use_fft = 1;
  80                                         break;
  81                                 case 'N':
  82                                         use_fft = 0;
  83                                         break;
  84                                 case 'e':
  85                                         fftscore = 0;
  86                                         break;
  87                                 case 'Q':
  88                                         alg = 'Q';
  89                                         break;
  90                                 case 'H':
  91                                         alg = 'H';
  92                                         break;
  93                                 case 'A':
  94                                         alg = 'A';
  95                                         break;
  96                                 case 'M':
  97                                         alg = 'M';
  98                                         break;
  99                                 case 'd':
 100                                         disp = 1;
 101                                         break;
 102                                 case 'O':
 103                                         outgap = 0;
 104                                         break;
 105                                 case 'a':
 106                                         fmodel = 1;
 107                                         break;
 108                                 case 'u':
 109                                         tbrweight = 0;
 110                                         break;
 111                                 case 'z':
 112                                         fftThreshold = atoi( *++argv );
 113                                         --argc;
 114                                         goto nextoption;
 115                                 case 'w':
 116                                         fftWinSize = atoi( *++argv );
 117                                         --argc;
 118                                         goto nextoption;
 119                                 case 'Z':
 120                                         checkC = 1;
 121                                         break;
 122                 case 'f':
 123                     ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 124                     fprintf( stderr, "ppenalty = %d\n", ppenalty );
 125                     --argc;
 126                     goto nextoption;
 127                 case 'g':
 128                     ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
 129                     fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
 130                     --argc;
 131                     goto nextoption;
 132                 case 'h':
 133                     poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
 134                     fprintf( stderr, "poffset = %d\n", poffset );
 135                     --argc;
 136                     goto nextoption;
 137                 case 'k':
 138                     kimuraR = atoi( *++argv );
 139                     fprintf( stderr, "kappa = %d\n", kimuraR );
 140                     --argc;
 141                     goto nextoption;
 142                                 case 'b':
 143                                         nblosum = atoi( *++argv );
 144                                         scoremtx = 1;
 145                                         fprintf( stderr, "blosum %d\n", nblosum );
 146                                         --argc;
 147                                         goto nextoption;
 148                                 case 'j':
 149                                         pamN = atoi( *++argv );
 150                                         scoremtx = 0;
 151                                         TMorJTT = JTT;
 152                                         fprintf( stderr, "jtt %d\n", pamN );
 153                                         --argc;
 154                                         goto nextoption;
 155                                 case 'm':
 156                                         pamN = atoi( *++argv );
 157                                         scoremtx = 0;
 158                                         TMorJTT = TM;
 159                                         fprintf( stderr, "tm %d\n", pamN );
 160                                         --argc;
 161                                         goto nextoption;
 162                                 default:
 163                                         fprintf( stderr, "illegal option %c\n", c );
 164                                         argc = 0;
 165                                         break;
 166             }
 167                 }
 168                 nextoption:
 169                         ;
 170         }
 171     if( argc != 2 )
 172     {
 173         fprintf( stderr, "options: Check source file ! %c ?\n", c );
 174         exit( 1 );
 175     }
 176         fprintf( stderr, "tbitr = %d, tbrweight = %d, tbweight = %d\n", tbitr, tbrweight, tbweight );
 177 //      readOtherOptions( &ppid, &fftThreshold, &fftWinSize );
 178         return( argv );
 179
 180 }
 181
 182 void GroupAlign( int nseq1, int nseq2, char **name, int *nlen, char **seq, char **aseq, char **mseq1, char **mseq2, int ***topol, double **len, double *eff, int alloclen )
 183 {
 184         int i;
 185         int clus1, clus2;
 186         int s1, s2;
 187         float pscore;
 188         static char **name1, **name2;
 189         double *effarr = eff;
 190         double *effarr1 = NULL;
 191         double *effarr2 = NULL;
 192         static char *indication1, *indication2;
 193         float dumfl = 0.0;
 194         int intdum;
 195 #if DEBUG
 196         double time1, time2;
 197 #endif
 198
 199
 200 //      fprintf( stderr, "in GroupAlign fftWinSize   = %d\n", fftWinSize );
 201 //      fprintf( stderr, "in GroupAlign fftThreshold = %d\n", fftThreshold );
 202
 203         if( effarr1 == NULL )
 204         {
 205                 name1 = AllocateCharMtx( nseq1, B );
 206                 name2 = AllocateCharMtx( nseq2, B );
 207                 indication1 = AllocateCharVec( 150 );
 208                 indication2 = AllocateCharVec( 150 );
 209                 effarr1 = AllocateDoubleVec( njob );
 210                 effarr2 = AllocateDoubleVec( njob );
 211 #if 0
 212 #else
 213 #endif
 214         }
 215
 216         for( i=0; i<njob; i++ ) strcpy( aseq[i], seq[i] );
 217
 218
 219
 220         s1 = 0;
 221         s2 = nseq1;
 222 //      fprintf( stdout, "nseq1 = %d\n", nseq1 );
 223
 224
 225
 226         clus1 = conjuctionforgaln( 0, nseq1, aseq, mseq1, effarr1, effarr, name, name1, indication1 );
 227         clus2 = conjuctionforgaln( nseq1, njob, aseq, mseq2, effarr2, effarr, name, name2, indication2 );
 228 /*
 229         fprintf( stderr, "before align all\n" );
 230         display( aseq, njob );
 231         fprintf( stderr, "\n" );
 232         fprintf( stderr, "before align 1 %s \n", indication1 );
 233         display( mseq1, clus1 );
 234         fprintf( stderr, "\n" );
 235         fprintf( stderr, "before align 2 %s \n", indication2 );
 236         display( mseq2, clus2 );
 237         fprintf( stderr, "\n" );
 238 */
 239
 240         if( alg == 'Q' )
 241                 fprintf( stderr, "before realignment, score = %f\n", naiveQpairscore( nseq1, nseq2, mseq1, mseq2, effarr1, effarr2, penalty ) );
 242         else if( alg == 'H' )
 243                 fprintf( stderr, "before realignment, score = %f\n", naivepairscore( nseq1, nseq2, mseq1, mseq2, effarr1, effarr2, penalty ) );
 244
 245         commongappick( nseq1, mseq1 );
 246         commongappick( nseq2, mseq2 );
 247
 248
 249 #if DEBUG
 250         time1 = getrusage_sec();
 251         fprintf( stderr, "entering Falign\n" );
 252 #endif
 253         if( use_fft )
 254         {
 255                 if( alg == 'M' )
 256                         pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum );
 257                 else
 258                         pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
 259         }
 260         else
 261         {
 262                 if( alg == 'M' )
 263                         pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
 264                 else if( alg == 'H' )
 265                         pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
 266                 else if( alg == 'Q' )
 267                         pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
 268                 else
 269                         pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
 270         }
 271 #if DEBUG
 272                 time2 = getrusage_sec();
 273                 fprintf( stdout, "### %d - %d,  %f\n", clus1, clus2, time2-time1 );
 274                 fflush( stdout );
 275 #endif
 276
 277         if( alg == 'H' )
 278                 fprintf( stderr, "before realignment, score = %f\n", naivepairscore( nseq1, nseq2, mseq1, mseq2, effarr1, effarr2, penalty ) );
 279         else if( alg == 'Q' )
 280                 fprintf( stderr, "after realignment, score = %f\n", naiveQpairscore( nseq1, nseq2, mseq1, mseq2, effarr1, effarr2, penalty ) );
 281
 282
 283 /*
 284         fprintf( stderr, "after align 1 %s \n", indication1 );
 285         display( mseq1, clus1 );
 286         fprintf( stderr, "\n" );
 287         fprintf( stderr, "after align 2 %s \n", indication2 );
 288         display( mseq2, clus2 );
 289         fprintf( stderr, "\n" );
 290 */
 291
 292         fprintf( stderr, "group-to-group %s /%s     %f\n", indication1, indication2, pscore );
 293         if( disp ) display( aseq, njob );
 294         fprintf( stderr, "\n" );
 295
 296 /*
 297         trap = fopen( "pre", "r+" );
 298         if( !trap ) ErrorExit( 1 );
 299         WriteGapFill( trap, njob, name, nlen, aseq );
 300         fclose( trap );
 301         fprintf( stdout, "nseq1 = %d\n", nseq1 );
 302 */
 303 }
 304
 305 static void WriteOptions( FILE *fp )
 306 {
 307         fprintf( fp, "tree-base method\n" );
 308         if( tbweight == 0 ) fprintf( fp, "unweighted\n" );
 309         else if( tbweight == 3 ) fprintf( fp, "reversely weighted\n" );
 310         if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
 311         else if( scoremtx ==  1 ) fprintf( fp, "Dayhoff( machigai ga aru )\n" );
 312         else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
 313         else if( scoremtx == -1 ) fprintf( fp, "DNA\n" );
 314
 315         if( scoremtx )
 316                 fprintf( fp, "Gap Penalty = %d, %d\n", penalty, offset );
 317         else
 318                 fprintf( fp, "Gap Penalty = %d\n", penalty );
 319 }
 320
 321
 322 int main( int argc, char *argv[] )
 323 {
 324         char **argv2;
 325         static int  *nlen;
 326         static char **name, **seq;
 327         static char **seq1, **seq2;
 328         static char **mseq1, **mseq2;
 329         static char **aseq;
 330         static char **bseq;
 331         static double **pscore;
 332         static double *eff;
 333         int i, j, len1, len2;
 334         static int ***topol;
 335         static double **len;
 336         FILE *gp1, *gp2;
 337         char c;
 338         int nlenmax1, nlenmax2, nseq1, nseq2;
 339         int alloclen;
 340
 341         argv2 = arguments( argc, argv );
 342
 343         fprintf( stderr, "####### in galn\n" );
 344
 345         initFiles();
 346
 347         fprintf( stderr, "file1 = %s\n", argv2[0] );
 348         fprintf( stderr, "file2 = %s\n", argv2[1] );
 349
 350         gp1 = fopen( argv2[0], "r" ); if( !gp1 ) ErrorExit( "cannot open file1" );
 351         gp2 = fopen( argv2[1], "r" ); if( !gp2 ) ErrorExit( "cannot open file2" );
 352
 353 #if 0
 354         PreRead( gp1, &nseq1, &nlenmax1 );
 355         PreRead( gp2, &nseq2, &nlenmax2 );
 356 #else
 357     getnumlen( gp1 );
 358         nseq1 = njob; nlenmax1 = nlenmax;
 359     getnumlen( gp2 );
 360         nseq2 = njob; nlenmax2 = nlenmax;
 361 #endif
 362
 363         njob = nseq1 + nseq2;
 364         nlenmax = MAX( nlenmax1, nlenmax2 );
 365
 366         rewind( gp1 );
 367         rewind( gp2 );
 368
 369
 370         name = AllocateCharMtx( njob, B );
 371         nlen = AllocateIntVec( njob );
 372         seq1 = AllocateCharMtx( nseq1, nlenmax*3 );
 373         seq2 = AllocateCharMtx( nseq2, nlenmax*3 );
 374         seq  = AllocateCharMtx( njob, 1 );
 375         aseq = AllocateCharMtx( njob, nlenmax*3 );
 376         bseq = AllocateCharMtx( njob, nlenmax*3 );
 377         mseq1 = AllocateCharMtx( njob, 1 );
 378         mseq2 = AllocateCharMtx( njob, 1 );
 379         alloclen = nlenmax * 3;
 380
 381         topol = AllocateIntCub( njob, 2, njob );
 382         len = AllocateDoubleMtx( njob, 2 );
 383         pscore = AllocateDoubleMtx( njob, njob );
 384         eff = AllocateDoubleVec( njob );
 385
 386 #if 0
 387     njob=nseq2; FRead( gp2, name+nseq1, nlen+nseq1, seq2 );
 388         njob=nseq1; FRead( gp1, name, nlen, seq1 );
 389 #else
 390     njob=nseq2; readDataforgaln( gp2, name+nseq1, nlen+nseq1, seq2 );
 391         njob=nseq1; readDataforgaln( gp1, name, nlen, seq1 );
 392 #endif
 393         njob = nseq1 + nseq2;
 394
 395
 396 #if 0  // CHUUI
 397         commongappick( nseq1, seq1 );
 398         commongappick( nseq2, seq2 );
 399 #endif
 400
 401         for( i=0; i<nseq1; i++ ) seq[i] = seq1[i];
 402         for( i=nseq1; i<njob; i++ ) seq[i] = seq2[i-nseq1];
 403 /*
 404         Write( stdout, njob, name, nlen, seq );
 405 */
 406
 407     constants( njob, seq );
 408
 409     WriteOptions( trap_g );
 410
 411     c = seqcheck( seq );
 412     if( c )
 413     {
 414         fprintf( stderr, "Illeagal character %c\n", c );
 415         exit( 1 );
 416     }
 417     for( i=1; i<nseq1; i++ )
 418     {
 419         if( nlen[i] != nlen[0] )
 420             ErrorExit( "group1 is not aligned." );
 421     }
 422     for( i=nseq1+1;  i<njob; i++ )
 423     {
 424         if( nlen[i] != nlen[nseq1] )
 425             ErrorExit( "group2 is not aligned." );
 426     }
 427     if( tbutree == 0 )
 428         {
 429                 for( i=0; i<nseq1; i++ )
 430                 {
 431                         for( j=i+1; j<nseq1; j++ )
 432                         {
 433                                 pscore[i][j] = (double)substitution_hosei( seq[i], seq[j] );
 434 //                              fprintf( stderr, "%d-%d, %5.1f \n", i, j, pscore[i][j] );
 435                         }
 436                         for( j=nseq1; j<njob; j++ )
 437                         {
 438                                 pscore[i][j] = 3.0;
 439 //                              fprintf( stderr, "%d-%d, %5.1f \n", i, j, pscore[i][j] );
 440                         }
 441                 }
 442                 for( i=nseq1; i<njob-1; i++ )
 443                 {
 444                         for( j=i+1; j<njob; j++ )
 445                         {
 446                                 pscore[i][j] = (double)substitution_hosei( seq[i], seq[j] );
 447 //                              fprintf( stderr, "%d-%d, %5.1f \n", i, j, pscore[i][j] );
 448                         }
 449                 }
 450 //              fprintf( stderr, "\n" );
 451
 452
 453     }
 454         else
 455         {
 456                 fprintf( stderr, "Not supported\n" );
 457                 exit( 1 );
 458 #if 0
 459                 prep = fopen( "hat2", "r" );
 460                 if( prep == NULL ) ErrorExit( "Make hat2." );
 461                 readhat2( prep, njob, name, pscore );
 462                 fclose( prep );
 463 #endif
 464         }
 465         fprintf( stderr, "Constructing dendrogram ... " );
 466         if( treemethod == 'x' )
 467                 veryfastsupg( njob, pscore, topol, len );
 468         else
 469                 ErrorExit( "Incorrect tree\n" );
 470         fprintf( stderr, "done.\n" );
 471
 472         if( tbrweight )
 473         {
 474                 weight = 3;
 475                 counteff_simple( njob, topol, len, eff );
 476 //              for( i=0; i<njob; i++ ) fprintf( stderr, "eff[%d] = %f\n", i, eff[i] );
 477         }
 478         else
 479         {
 480                 for( i=0; i<njob; i++ ) eff[i] = 1.0;
 481         }
 482
 483         len1 = strlen( seq[0] );
 484         len2 = strlen( seq[nseq1] );
 485         if( len1 > 30000 || len2 > 30000 )
 486         {
 487                 fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
 488                 alg = 'M';
 489         }
 490
 491
 492
 493
 494         GroupAlign( nseq1, nseq2, name, nlen, seq, aseq, mseq1, mseq2, topol, len, eff, alloclen );
 495
 496 #if 0
 497         writePre( njob, name, nlen, aseq, 1 );
 498 #else
 499         writeDataforgaln( stdout, njob, name, nlen, aseq );
 500 #endif
 501
 502         SHOWVERSION;
 503         return( 0 );
 504 }