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