binaries/src/mafft/core/dndblast.c

   1 #include "mltaln.h"
   2 #include <sys/types.h>
   3 #include <unistd.h>
   4 #define DEBUG 0
   5 #define TEST 0
   6
   7
   8 int howmanyx( char *s )
   9 {
  10         int val = 0;
  11         if( scoremtx == -1 )
  12         {
  13                 do
  14                 {
  15                         if( !strchr( "atgcuATGCU", *s ) ) val++;
  16                 } while( *++s );
  17         }
  18         else
  19         {
  20                 do
  21                 {
  22                         if( !strchr( "ARNDCQEGHILKMFPSTWYV", *s ) ) val++;
  23                 } while( *++s );
  24         }
  25         return( val );
  26 }
  27
  28 void arguments( int argc, char *argv[] )
  29 {
  30         int c;
  31
  32         inputfile = NULL;
  33         disopt = 0;
  34         divpairscore = 0;
  35
  36     while( --argc > 0 && (*++argv)[0] == '-' )
  37         {
  38         while ( (c = *++argv[0]) )
  39                 {
  40             switch( c )
  41             {
  42                                 case 'i':
  43                                         inputfile = *++argv;
  44                                         fprintf( stderr, "inputfile = %s\n", inputfile );
  45                                         --argc;
  46                                         goto nextoption;
  47                                 case 'I':
  48                                         disopt = 1;
  49                                         break;
  50                 default:
  51                     fprintf( stderr, "illegal option %c\n", c );
  52                     argc = 0;
  53                     break;
  54             }
  55                 }
  56                 nextoption:
  57                         ;
  58
  59         }
  60     if( argc != 0 )
  61     {
  62         fprintf( stderr, "options: -i\n" );
  63         exit( 1 );
  64     }
  65 }
  66
  67 int main( int argc, char *argv[] )
  68 {
  69         int ktuple;
  70         int i, j;
  71         FILE *infp;
  72         FILE *hat2p;
  73         FILE *hat3p;
  74         char **seq = NULL; // by D.Mathog
  75         char **seq1;
  76         static char name[M][B];
  77         static char name1[M][B];
  78         static int nlen1[M];
  79         double **mtx;
  80         double **mtx2;
  81         static int nlen[M];
  82         char b[B];
  83         double max;
  84         char com[1000];
  85         int opt[M];
  86         int res;
  87         char *home;
  88         char queryfile[B];
  89         char datafile[B];
  90         char fastafile[B];
  91         char hat2file[B];
  92         int pid = (int)getpid();
  93         LocalHom **localhomtable, *tmpptr;
  94 #if 1
  95         home = getenv( "HOME" );
  96 #else /* $HOME wo tsukau to fasta ni watasu hikisuu ga afureru */
  97         home = NULL;
  98 #endif
  99
 100 #if DEBUG
 101         if( home ) fprintf( stderr, "home = %s\n", home );
 102 #endif
 103         if( !home ) home = "";
 104         sprintf( queryfile, "%s/tmp/query-%d", home, pid );
 105         sprintf( datafile, "%s/tmp/data-%d", home, pid );
 106         sprintf( fastafile, "%s/tmp/fasta-%d", home, pid );
 107         sprintf( hat2file, "hat2-%d", pid );
 108
 109
 110         arguments( argc, argv );
 111
 112         if( inputfile )
 113         {
 114                 infp = fopen( inputfile, "r" );
 115                 if( !infp )
 116                 {
 117                         fprintf( stderr, "Cannot open %s\n", inputfile );
 118                         exit( 1 );
 119                 }
 120         }
 121         else
 122                 infp = stdin;
 123 #if 0
 124         PreRead( infp, &njob, &nlenmax );
 125 #else
 126         dorp = NOTSPECIFIED;
 127         getnumlen( infp );
 128 #endif
 129
 130         if( dorp == 'd' )
 131         {
 132                 scoremtx = -1;
 133                 pamN = NOTSPECIFIED;
 134         }
 135         else
 136         {
 137                 nblosum = 62;
 138                 scoremtx = 1;
 139         }
 140         constants( njob, seq );
 141
 142         rewind( infp );
 143
 144         seq = AllocateCharMtx( njob, nlenmax+1 );
 145         seq1 = AllocateCharMtx( 2, nlenmax+1 );
 146         mtx = AllocateDoubleMtx( njob, njob );
 147         mtx2 = AllocateDoubleMtx( njob, njob );
 148         localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 149         for( i=0; i<njob; i++)
 150         {
 151                 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
 152                 for( j=0; j<njob; j++)
 153                 {
 154                         localhomtable[i][j].start1 = -1;
 155                         localhomtable[i][j].end1 = -1;
 156                         localhomtable[i][j].start2 = -1;
 157                         localhomtable[i][j].end2 = -1;
 158                         localhomtable[i][j].opt = -1.0;
 159                         localhomtable[i][j].next = NULL;
 160
 161                 }
 162     }
 163
 164 #if 0
 165         FRead( infp, name, nlen, seq );
 166 #else
 167         readData( infp, name, nlen, seq );
 168 #endif
 169         fclose( infp );
 170
 171         if( scoremtx == -1 ) ktuple = 6;
 172         else                 ktuple = 1;
 173
 174         for( i=0; i<njob; i++ )
 175         {
 176                 gappick0( seq1[0], seq[i] );
 177                 strcpy( seq[i], seq1[0] );
 178         }
 179         for( j=0; j<njob; j++ )
 180         {
 181                 sprintf( name1[j], "+==========+%d                      ", j );
 182                 nlen1[j] = nlen[j];
 183         }
 184
 185         for( i=0; i<njob; i++ )
 186         {
 187 //              fprintf( stderr, "###  i = %d\n", i );
 188
 189                 if( i % 10 == 0 )
 190                 {
 191                         fprintf( stderr, "formatting .. " );
 192                         hat2p = fopen( datafile, "w" );
 193                         if( !hat2p ) ErrorExit( "Cannot open datafile." );
 194                         WriteForFasta( hat2p, njob-i, name1+i, nlen1+i, seq+i );
 195                         fclose( hat2p );
 196
 197                         if( scoremtx == -1 )
 198                                 sprintf( com, "formatdb  -p f -i %s -o F", datafile );
 199                         else
 200                                 sprintf( com, "formatdb  -i %s -o F", datafile );
 201                         system( com );
 202                         fprintf( stderr, "done.\n" );
 203                 }
 204
 205                 hat2p = fopen( queryfile, "w" );
 206                 if( !hat2p ) ErrorExit( "Cannot open queryfile." );
 207                 WriteForFasta( hat2p, 1, name1+i, nlen+i, seq+i );
 208                 fclose( hat2p );
 209
 210
 211                 if( scoremtx == -1 )
 212                         sprintf( com, "blastall -e 1e10 -p blastn -m 7  -i %s -d %s >  %s", queryfile, datafile, fastafile );
 213                 else
 214                         sprintf( com, "blastall -G 10 -E 1 -e 1e10 -p blastp -m 7  -i %s -d %s >  %s", queryfile, datafile, fastafile );
 215                 res = system( com );
 216                 if( res ) ErrorExit( "error in fasta" );
 217
 218
 219                 hat2p = fopen( fastafile, "r" );
 220                 if( hat2p == NULL )
 221                         ErrorExit( "file 'fasta.$$' does not exist\n" );
 222                 res = ReadBlastm7( hat2p, mtx[i], i, name1, localhomtable[i] );
 223                 fclose( hat2p );
 224
 225 #if 0
 226                 for( j=0; j<njob; j++ )
 227                 {
 228                         for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
 229                         {
 230                                 if( tmpptr->opt == -1.0 ) continue;
 231 //                              fprintf( stderr, "%d %d %d %6.3f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->next );
 232                         }
 233                 }
 234 #endif
 235
 236                 if( res < njob-i+i%10 )
 237                 {
 238                         fprintf( stderr, "WARNING: count (blast) = %d < %d\n", res, njob-i+i%10 );
 239                 }
 240
 241
 242 #if 0
 243                 {
 244                         int ii, jj;
 245                         if( i < njob-1 ) for( jj=i; jj<i+5; jj++ )
 246                                 fprintf( stdout, "mtx[%d][%d] = %f\n", i+1, jj+1, mtx[i][jj] );
 247                 }
 248 #endif
 249                 fprintf( stderr, "query : %4d / %d\n", i+1, njob );
 250         }
 251
 252 #if 1
 253         fprintf( stderr, "##### writing hat3\n" );
 254         hat3p = fopen( "hat3", "w" );
 255         if( !hat3p ) ErrorExit( "Cannot open hat3." );
 256         for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ )
 257         {
 258 //              fprintf( stderr, "mtx[%d][%d] = %f, mtx[%d][%d] = %f\n", i, j, mtx[i][j], j, i, mtx[j][i] );
 259                 if( i == j ) continue;
 260                 if( mtx[i][j] == mtx[j][i] && i > j ) continue;
 261                 if( mtx[j][i] > mtx[i][j] ) continue;
 262                 for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
 263                 {
 264                         if( tmpptr->opt == -1.0 ) continue;
 265                         fprintf( hat3p, "%d %d %d %6.3f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next );
 266                 }
 267         }
 268         fclose( hat3p );
 269 #endif
 270
 271         for( i=0; i<njob; i++ )
 272         {
 273 //              fprintf( stderr, "###  i = %d\n", i );
 274                 hat2p = fopen( datafile, "w" );
 275                 if( !hat2p ) ErrorExit( "Cannot open datafile." );
 276                 WriteForFasta( hat2p, njob-i, name1+i, nlen1+i, seq+i );
 277                 fclose( hat2p );
 278
 279 //              seq1[0] = seq[i];
 280 //              nlen1[0] = nlen[i];
 281
 282                 hat2p = fopen( queryfile, "w" );
 283                 if( !hat2p ) ErrorExit( "Cannot open queryfile." );
 284                 WriteForFasta( hat2p, 1, name1+i, nlen+i, seq+i );
 285                 fclose( hat2p );
 286
 287
 288                 if( scoremtx == -1 )
 289                         sprintf( com, "fasta34 -z3 -m10  -n -Q  -b%d -E%d -d%d %s %s %d > %s", M, M, 0, queryfile, datafile, ktuple, fastafile );
 290                 else
 291                         sprintf( com, "fasta34 -z3 -m10  -Q  -b%d -E%d -d%d %s %s %d > %s", M, M, 0, queryfile, datafile, ktuple, fastafile );
 292                 res = system( com );
 293                 if( res ) ErrorExit( "error in fasta" );
 294
 295
 296                 hat2p = fopen( fastafile, "r" );
 297                 if( hat2p == NULL )
 298                         ErrorExit( "file 'fasta.$$' does not exist\n" );
 299                 res = ReadFasta34noalign( hat2p, mtx[i], i, name1, localhomtable[i] );
 300                 fclose( hat2p );
 301                 if( res < njob - i )
 302                 {
 303                         fprintf( stderr, "count (fasta34 -z 3) = %d\n", res );
 304                         exit( 1 );
 305                 }
 306
 307
 308                 if( i == 0 )
 309                         for( j=0; j<njob; j++ ) opt[j] = (int)mtx[0][j];
 310
 311
 312 #if 0
 313                 {
 314                         int ii, jj;
 315                         if( i < njob-1 ) for( jj=i; jj<i+5; jj++ )
 316                                 fprintf( stdout, "mtx[%d][%d] = %f\n", i+1, jj+1, mtx[i][jj] );
 317                 }
 318 #endif
 319                 fprintf( stderr, "query : %4d\r", i+1 );
 320         }
 321
 322
 323
 324
 325         for( i=0; i<njob; i++ )
 326         {
 327                 max = mtx[i][i];
 328                 if( max == 0.0 )
 329                 {
 330                         for( j=0; j<njob; j++ )
 331                                 mtx2[i][j] = 2.0;
 332                 }
 333                 else
 334                 {
 335                         for( j=0; j<njob; j++ )
 336                         {
 337                                 mtx2[i][j] = ( max - mtx[MIN(i,j)][MAX(i,j)] ) / max * 2.0;
 338 //                              fprintf( stdout, "max = %f, mtx[%d][%d] = %f -> %f\n", max, i+1, j+1, mtx[i][j], mtx2[i][j] );
 339                         }
 340                 }
 341         }
 342         for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 343         {
 344 //              fprintf( stdout, "mtx2[%d][%d] = %f, %f\n", i+1, j+1, mtx2[i][j], mtx2[j][i] );
 345                 mtx2[i][j] = MIN( mtx2[i][j], mtx2[j][i] );
 346         }
 347
 348 #if 0
 349         {
 350                 int ii, jj;
 351                 if( i < njob-1 ) for( jj=i+1; jj<njob; jj++ )
 352                         fprintf( stderr, "mtx2[][] = %f\n", mtx2[i][jj] );
 353         }
 354 #endif
 355
 356         for( i=0; i<njob; i++ ) name[i][0] = '=';
 357
 358         if( disopt )
 359         {
 360                 strcpy( b, name[0] );
 361                 sprintf( name[0], "=query====lgth=%04d-%04d %.*s", nlen[0], howmanyx( seq[0] ), B-30, b );
 362 #if 0
 363                 strins(  b, name[0] );
 364 #endif
 365                 for( i=1; i<njob; i++ )
 366                 {
 367                         strcpy( b, name[i] );
 368                         sprintf( name[i], "=opt=%04d=lgth=%04d-%04d %.*s", opt[i], nlen[i], howmanyx( seq[i] ), B-30, b );
 369 #if 0
 370                         strins( b, name[i] );
 371 #endif
 372                 }
 373         }
 374
 375         hat2p = fopen( hat2file, "w" );
 376         if( !hat2p ) ErrorExit( "Cannot open hat2." );
 377         WriteHat2( hat2p, njob, name, mtx2 );
 378         fclose( hat2p );
 379
 380
 381         sprintf( com, "/bin/rm %s %s %s", queryfile, datafile, fastafile );
 382         system( com );
 383
 384 #if 0
 385         sprintf( com, ALNDIR "/supgsdl < %s", hat2file );
 386         res = system( com );
 387         if( res ) ErrorExit( "error in spgsdl" );
 388 #endif
 389
 390         sprintf( com, "mv %s hat2", hat2file );
 391         res = system( com );
 392         if( res ) ErrorExit( "error in mv" );
 393
 394         SHOWVERSION;
 395         exit( 0 );
 396 }