binaries/src/mafft/core/galn.c

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