binaries/src/mafft/core/dndfast4.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         disopt = 0;
  33
  34     while( --argc > 0 && (*++argv)[0] == '-' )
  35         while ( c = *++argv[0] )
  36             switch( c )
  37             {
  38                                 case 'i':
  39                                         disopt = 1;
  40                                         break;
  41                 default:
  42                     fprintf( stderr, "illegal option %c\n", c );
  43                     argc = 0;
  44                     break;
  45             }
  46     if( argc != 0 )
  47     {
  48         fprintf( stderr, "options: -i\n" );
  49         exit( 1 );
  50     }
  51 }
  52
  53 int main( int argc, char *argv[] )
  54 {
  55         int ktuple;
  56         int i, j;
  57         FILE *hat2p;
  58         char **seq;
  59         char **seq1;
  60         static char name[M][B];
  61         static char name1[M][B];
  62         static int nlen1[M];
  63         double **mtx;
  64         double **mtx2;
  65         static int nlen[M];
  66         char b[B];
  67         double max;
  68         char com[B];
  69         int opt[M];
  70         int res;
  71         char *home;
  72         char queryfile[B];
  73         char datafile[B];
  74         char fastafile[B];
  75         char hat2file[B];
  76         int pid = (int)getpid();
  77 #if 0
  78         home = getenv( "HOME" );
  79 #else /* $HOME wo tsukau to fasta ni watasu hikisuu ga afureru */
  80         home = NULL;
  81 #endif
  82
  83 #if DEBUG
  84         if( home ) fprintf( stderr, "home = %s\n", home );
  85 #endif
  86         if( !home ) home = "";
  87         sprintf( queryfile, "%s/tmp/query-%d\0", home, pid );
  88         sprintf( datafile, "%s/tmp/data-%d\0", home, pid );
  89         sprintf( fastafile, "%s/tmp/fasta-%d\0", home, pid );
  90         sprintf( hat2file, "hat2-%d\0", pid );
  91
  92         arguments( argc, argv );
  93 #if 0
  94         PreRead( stdin, &njob, &nlenmax );
  95 #else
  96         getnumlen( stdin );
  97 #endif
  98         rewind( stdin );
  99
 100         seq = AllocateCharMtx( njob, nlenmax+1 );
 101         seq1 = AllocateCharMtx( 2, nlenmax+1 );
 102         mtx = AllocateDoubleMtx( njob, njob );
 103         mtx2 = AllocateDoubleMtx( njob, njob );
 104
 105 #if 0
 106         FRead( stdin, name, nlen, seq );
 107 #else
 108         readData( stdin, name, nlen, seq );
 109 #endif
 110         if( scoremtx == -1 ) ktuple = 6;
 111         else                 ktuple = 1;
 112
 113         for( i=0; i<njob; i++ )
 114         {
 115                 gappick0( seq1[0], seq[i] );
 116                 strcpy( seq[i], seq1[0] );
 117         }
 118         for( j=0; j<njob; j++ )
 119         {
 120                 sprintf( name1[j], "+==========+%d                      \0", j );
 121                 nlen1[j] = nlen[j];
 122         }
 123         hat2p = fopen( datafile, "w" );
 124         if( !hat2p ) ErrorExit( "Cannot open datafile." );
 125         WriteForFasta( hat2p, njob, name1, nlen1, seq );
 126         fclose( hat2p );
 127
 128         for( i=0; i<njob; i++ )
 129         {
 130
 131                 hat2p = fopen( datafile, "w" );
 132                 if( !hat2p ) ErrorExit( "Cannot open datafile." );
 133                 WriteForFasta( hat2p, njob-i, name1+i, nlen1+i, seq+i );
 134                 fclose( hat2p );
 135
 136
 137                 seq1[0] = seq[i];
 138                 nlen1[0] = nlen[i];
 139
 140                 hat2p = fopen( queryfile, "w" );
 141                 if( !hat2p ) ErrorExit( "Cannot open queryfile." );
 142                 WriteForFasta( hat2p, 1, name1+i, nlen1, seq1 );
 143                 fclose( hat2p );
 144
 145                 if( scoremtx == -1 )
 146                         sprintf( com, "fasta3 -n -Q -h -b%d -E%d -d%d %s %s %d > %s\0", M, M, 0, queryfile, datafile, ktuple, fastafile );
 147                 else
 148                         sprintf( com, "fasta3 -Q -h -b%d -E%d -d%d %s %s %d > %s\0", M, M, 0, queryfile, datafile, ktuple, fastafile );
 149                 res = system( com );
 150                 if( res ) ErrorExit( "error in fasta" );
 151
 152                 hat2p = fopen( fastafile, "r" );
 153                 if( hat2p == NULL )
 154                         ErrorExit( "file 'fasta.$$' does not exist\n" );
 155                 ReadFasta3( hat2p, mtx[i], njob-i, name1 );
 156
 157                 if( i == 0 )
 158                         for( j=0; j<njob; j++ ) opt[j] = (int)mtx[0][j];
 159
 160                 fclose( hat2p );
 161
 162 #if 1
 163                 {
 164                         int ii, jj;
 165                         if( i < njob-1 ) for( jj=i; jj<i+5; jj++ )
 166                                 fprintf( stdout, "mtx[%d][%d] = %f\n", i+1, jj+1, mtx[i][jj] );
 167                 }
 168 #endif
 169                 fprintf( stderr, "query : %#4d\n", i+1 );
 170         }
 171
 172         for( i=0; i<njob; i++ )
 173         {
 174                 max = mtx[i][i];
 175                 if( max == 0.0 )
 176                 {
 177                         for( j=0; j<njob; j++ )
 178                                 mtx2[i][j] = 2.0;
 179                 }
 180                 else
 181                 {
 182                         for( j=0; j<njob; j++ )
 183                         {
 184                                 mtx2[i][j] = ( max - mtx[MIN(i,j)][MAX(i,j)] ) / max * 2.0;
 185 //                              fprintf( stdout, "max = %f, mtx[%d][%d] = %f -> %f\n", max, i+1, j+1, mtx[i][j], mtx2[i][j] );
 186                         }
 187                 }
 188         }
 189         for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 190         {
 191 //              fprintf( stdout, "mtx2[%d][%d] = %f, %f\n", i+1, j+1, mtx2[i][j], mtx2[j][i] );
 192                 mtx2[i][j] = MIN( mtx2[i][j], mtx2[j][i] );
 193         }
 194 #if 0
 195         {
 196                 int ii, jj;
 197                 if( i < njob-1 ) for( jj=i+1; jj<njob; jj++ )
 198                         fprintf( stderr, "mtx2[][] = %f\n", mtx2[i][jj] );
 199         }
 200 #endif
 201
 202         for( i=0; i<njob; i++ ) name[i][0] = '=';
 203
 204         if( disopt )
 205         {
 206                 strcpy( b, name[0] );
 207                 sprintf( name[0], "=query====lgth=%#04d-%04d %.*s\0", nlen[0], howmanyx( seq[0] ), B-30, b );
 208 #if 0
 209                 strins(  b, name[0] );
 210 #endif
 211                 for( i=1; i<njob; i++ )
 212                 {
 213                         strcpy( b, name[i] );
 214                         sprintf( name[i], "=opt=%#04d=lgth=%#04d-%04d %.*s\0", opt[i], nlen[i], howmanyx( seq[i] ), B-30, b );
 215 #if 0
 216                         strins( b, name[i] );
 217 #endif
 218                 }
 219         }
 220
 221         hat2p = fopen( hat2file, "w" );
 222         if( !hat2p ) ErrorExit( "Cannot open hat2." );
 223         WriteHat2( hat2p, njob, name, mtx2 );
 224         fclose( hat2p );
 225
 226         sprintf( com, "/bin/rm %s %s %s", queryfile, datafile, fastafile );
 227         system( com );
 228
 229 #if 0
 230         sprintf( com, ALNDIR "/supgsdl < %s\0", hat2file );
 231         res = system( com );
 232         if( res ) ErrorExit( "error in spgsdl" );
 233 #endif
 234
 235         sprintf( com, "mv %s hat2", hat2file );
 236         res = system( com );
 237         if( res ) ErrorExit( "error in mv" );
 238
 239         SHOWVERSION;
 240         exit( 0 );
 241 }