binaries/src/mafft/core/dndfast7.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         swopt = "";
  36
  37     while( --argc > 0 && (*++argv)[0] == '-' )
  38         {
  39         while ( (c = *++argv[0]) )
  40                 {
  41             switch( c )
  42             {
  43                                 case 'i':
  44                                         inputfile = *++argv;
  45                                         fprintf( stderr, "inputfile = %s\n", inputfile );
  46                                         --argc;
  47                                         goto nextoption;
  48                                 case 'I':
  49                                         disopt = 1;
  50                                         break;
  51                                 case 'A':
  52                                         swopt = "-A";
  53                                         break;
  54                 default:
  55                     fprintf( stderr, "illegal option %c\n", c );
  56                     argc = 0;
  57                     break;
  58             }
  59                 }
  60                 nextoption:
  61                         ;
  62         }
  63     if( argc != 0 )
  64     {
  65         fprintf( stderr, "options: -i\n" );
  66         exit( 1 );
  67     }
  68 }
  69
  70 int main( int argc, char *argv[] )
  71 {
  72         int ktuple;
  73         int i, j;
  74         FILE *hat2p;
  75         FILE *hat3p;
  76         FILE *infp;
  77         char **seq = NULL; // by D.Mathog
  78         char **seq1;
  79         static char name[M][B];
  80         static char name1[M][B];
  81         static int nlen1[M];
  82         double **mtx;
  83         double **mtx2;
  84         static int nlen[M];
  85         char b[B];
  86         double max;
  87         char com[1000];
  88         int opt[M];
  89         int res;
  90         char *home;
  91         char *fastapath;
  92         char queryfile[B];
  93         char datafile[B];
  94         char fastafile[B];
  95         char hat2file[B];
  96         int pid = (int)getpid();
  97         LocalHom **localhomtable, *tmpptr;
  98 #if 0
  99         home = getenv( "HOME" );
 100 #else /* $HOME wo tsukau to fasta ni watasu hikisuu ga afureru */
 101         home = NULL;
 102 #endif
 103         fastapath = getenv( "FASTA_4_MAFFT" );
 104         if( !fastapath )
 105                 fastapath = "fasta34";
 106
 107 #if DEBUG
 108         if( home ) fprintf( stderr, "home = %s\n", home );
 109 #endif
 110         if( !home ) home = "";
 111         sprintf( queryfile, "%s/tmp/query-%d", home, pid );
 112         sprintf( datafile, "%s/tmp/data-%d", home, pid );
 113         sprintf( fastafile, "%s/tmp/fasta-%d", home, pid );
 114         sprintf( hat2file, "hat2-%d", pid );
 115
 116
 117         arguments( argc, argv );
 118
 119         if( inputfile )
 120         {
 121                 infp = fopen( inputfile, "r" );
 122                 if( !infp )
 123                 {
 124                         fprintf( stderr, "Cannot open %s\n", inputfile );
 125                         exit( 1 );
 126                 }
 127         }
 128         else
 129                 infp = stdin;
 130
 131
 132
 133 #if 0
 134         PreRead( stdin, &njob, &nlenmax );
 135 #else
 136         dorp = NOTSPECIFIED;
 137         getnumlen( infp );
 138 #endif
 139
 140         if( dorp == 'd' )
 141         {
 142                 scoremtx = -1;
 143                 pamN = NOTSPECIFIED;
 144         }
 145         else
 146         {
 147                 nblosum = 62;
 148                 scoremtx = 1;
 149         }
 150         constants( njob, seq );
 151
 152         rewind( infp );
 153
 154         seq = AllocateCharMtx( njob, nlenmax+1 );
 155         seq1 = AllocateCharMtx( 2, nlenmax+1 );
 156         mtx = AllocateDoubleMtx( njob, njob );
 157         mtx2 = AllocateDoubleMtx( njob, njob );
 158         localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 159         for( i=0; i<njob; i++)
 160         {
 161                 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
 162                 for( j=0; j<njob; j++)
 163                 {
 164                         localhomtable[i][j].start1 = -1;
 165                         localhomtable[i][j].end1 = -1;
 166                         localhomtable[i][j].start2 = -1;
 167                         localhomtable[i][j].end2 = -1;
 168                         localhomtable[i][j].opt = -1.0;
 169                         localhomtable[i][j].next = NULL;
 170                 }
 171     }
 172
 173 #if 0
 174         FRead( infp, name, nlen, seq );
 175 #else
 176         readData( infp, name, nlen, seq );
 177 #endif
 178         fclose( infp );
 179
 180         if( scoremtx == -1 ) ktuple = 6;
 181         else                 ktuple = 1;
 182
 183         for( i=0; i<njob; i++ )
 184         {
 185                 gappick0( seq1[0], seq[i] );
 186                 strcpy( seq[i], seq1[0] );
 187         }
 188         for( j=0; j<njob; j++ )
 189         {
 190                 sprintf( name1[j], "+==========+%d                      ", j );
 191                 nlen1[j] = nlen[j];
 192         }
 193         hat2p = fopen( datafile, "w" );
 194         if( !hat2p ) ErrorExit( "Cannot open datafile." );
 195         WriteForFasta( hat2p, njob, name1, nlen1, seq );
 196         fclose( hat2p );
 197
 198         for( i=0; i<njob; i++ )
 199         {
 200 //              fprintf( stderr, "###  i = %d\n", i );
 201                 hat2p = fopen( datafile, "w" );
 202                 if( !hat2p ) ErrorExit( "Cannot open datafile." );
 203                 WriteForFasta( hat2p, njob-i, name1+i, nlen1+i, seq+i );
 204                 fclose( hat2p );
 205
 206                 seq1[0] = seq[i];
 207                 nlen1[0] = nlen[i];
 208
 209                 hat2p = fopen( queryfile, "w" );
 210                 if( !hat2p ) ErrorExit( "Cannot open queryfile." );
 211                 WriteForFasta( hat2p, 1, name1+i, nlen1, seq1 );
 212                 fclose( hat2p );
 213
 214
 215                 if( scoremtx == -1 )
 216                         sprintf( com, "%s %s -z3 -m10  -n -Q  -b%d -E%d -d%d %s %s %d > %s", fastapath, swopt, M, M, M, queryfile, datafile, ktuple, fastafile );
 217                 else
 218                         sprintf( com, "%s %s -z3 -m10  -Q  -b%d -E%d -d%d %s %s %d > %s", fastapath, swopt, M, M, M, queryfile, datafile, ktuple, fastafile );
 219                 res = system( com );
 220                 if( res ) ErrorExit( "error in fasta" );
 221
 222
 223
 224                 hat2p = fopen( fastafile, "r" );
 225                 if( hat2p == NULL )
 226                         ErrorExit( "file 'fasta.$$' does not exist\n" );
 227                 if( scoremtx == -1 )
 228                         res = ReadFasta34m10_nuc( hat2p, mtx[i], i, name1, localhomtable[i] );
 229                 else
 230                         res = ReadFasta34m10( hat2p, mtx[i], i, name1, localhomtable[i] );
 231                 fclose( hat2p );
 232
 233                 if( res < njob - i )
 234                 {
 235                         fprintf( stderr, "count (fasta34 -z 3) = %d\n", res );
 236                         exit( 1 );
 237                 }
 238
 239
 240                 if( i == 0 )
 241                         for( j=0; j<njob; j++ ) opt[j] = (int)mtx[0][j];
 242
 243
 244 #if 0
 245                 {
 246                         int ii, jj;
 247                         if( i < njob-1 ) for( jj=i; jj<i+5; jj++ )
 248                                 fprintf( stdout, "mtx[%d][%d] = %f\n", i+1, jj+1, mtx[i][jj] );
 249                 }
 250 #endif
 251                 fprintf( stderr, "query : %4d / %5d\r", i+1, njob );
 252         }
 253
 254         for( i=0; i<njob; i++ )
 255         {
 256                 max = mtx[i][i];
 257                 if( max == 0.0 )
 258                 {
 259                         for( j=0; j<njob; j++ )
 260                                 mtx2[i][j] = 2.0;
 261                 }
 262                 else
 263                 {
 264                         for( j=0; j<njob; j++ )
 265                         {
 266 //                              fprintf( stderr, "##### mtx[%d][%d] = %f\n", i, j, mtx[i][j] );
 267                                 mtx2[i][j] = ( max - mtx[MIN(i,j)][MAX(i,j)] ) / max * 2.0;
 268 //                              fprintf( stdout, "max = %f, mtx[%d][%d] = %f -> %f\n", max, i+1, j+1, mtx[i][j], mtx2[i][j] );
 269                         }
 270                 }
 271         }
 272         for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 273         {
 274 //              fprintf( stdout, "mtx2[%d][%d] = %f, %f\n", i+1, j+1, mtx2[i][j], mtx2[j][i] );
 275                 mtx2[i][j] = MIN( mtx2[i][j], mtx2[j][i] );
 276         }
 277
 278 #if 0
 279         {
 280                 int ii, jj;
 281                 if( i < njob-1 ) for( jj=i+1; jj<njob; jj++ )
 282                         fprintf( stderr, "mtx2[][] = %f\n", mtx2[i][jj] );
 283         }
 284 #endif
 285
 286         for( i=0; i<njob; i++ ) name[i][0] = '=';
 287
 288         if( disopt )
 289         {
 290                 strcpy( b, name[0] );
 291                 sprintf( name[0], "=query====lgth=%04d-%04d %.*s", nlen[0], howmanyx( seq[0] ), B-30, b );
 292 #if 0
 293                 strins(  b, name[0] );
 294 #endif
 295                 for( i=1; i<njob; i++ )
 296                 {
 297                         strcpy( b, name[i] );
 298                         sprintf( name[i], "=opt=%04d=lgth=%04d-%04d %.*s", opt[i], nlen[i], howmanyx( seq[i] ), B-30, b );
 299 #if 0
 300                         strins( b, name[i] );
 301 #endif
 302                 }
 303         }
 304
 305         hat2p = fopen( hat2file, "w" );
 306         if( !hat2p ) ErrorExit( "Cannot open hat2." );
 307         WriteHat2( hat2p, njob, name, mtx2 );
 308         fclose( hat2p );
 309
 310 #if 1
 311         fprintf( stderr, "##### writing hat3\n" );
 312         hat3p = fopen( "hat3", "w" );
 313         if( !hat3p ) ErrorExit( "Cannot open hat3." );
 314         for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 315         {
 316                 for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
 317                 {
 318                         if( tmpptr->opt == -1.0 ) continue;
 319                         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 );
 320                 }
 321         }
 322         fclose( hat3p );
 323 #endif
 324
 325         sprintf( com, "/bin/rm %s %s %s", queryfile, datafile, fastafile );
 326         system( com );
 327
 328 #if 0
 329         sprintf( com, ALNDIR "/supgsdl < %s", hat2file );
 330         res = system( com );
 331         if( res ) ErrorExit( "error in spgsdl" );
 332 #endif
 333
 334         sprintf( com, "mv %s hat2", hat2file );
 335         res = system( com );
 336         if( res ) ErrorExit( "error in mv" );
 337
 338         SHOWVERSION;
 339         exit( 0 );
 340 }