binaries/src/mafft/core/dndpre2.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 1
   6
   7
   8 void arguments( int argc, char *argv[] )
   9 {
  10     int c;
  11
  12         nblosum = 62;
  13         fmodel = 0;
  14         calledByXced = 0;
  15         scoremtx = 0;
  16         dorp = NOTSPECIFIED;
  17         ppenalty = NOTSPECIFIED;
  18         ppenalty_ex = NOTSPECIFIED;
  19         poffset = NOTSPECIFIED;
  20         kimuraR = NOTSPECIFIED;
  21         pamN = NOTSPECIFIED;
  22         geta2 = GETA2;
  23         fftWinSize = NOTSPECIFIED;
  24         fftThreshold = NOTSPECIFIED;
  25
  26     while( --argc > 0 && (*++argv)[0] == '-' )
  27         {
  28         while ( (c = *++argv[0]) )
  29                 {
  30             switch( c )
  31             {
  32                                 case 'f':
  33                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
  34                                         fprintf( stderr, "ppenalty = %d\n", ppenalty );
  35                                         --argc;
  36                                         goto nextoption;
  37                                 case 'g':
  38                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
  39                                         fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
  40                                         --argc;
  41                                         goto nextoption;
  42                                 case 'h':
  43                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
  44                                         fprintf( stderr, "poffset = %d\n", poffset );
  45                                         --argc;
  46                                         goto nextoption;
  47                                 case 'k':
  48                                         kimuraR = atoi( *++argv );
  49                                         fprintf( stderr, "kimuraR = %d\n", kimuraR );
  50                                         --argc;
  51                                         goto nextoption;
  52                                 case 'b':
  53                                         nblosum = atoi( *++argv );
  54                                         scoremtx = 1;
  55                                         fprintf( stderr, "blosum %d\n", nblosum );
  56                                         --argc;
  57                                         goto nextoption;
  58                                 case 'j':
  59                                         pamN = atoi( *++argv );
  60                                         scoremtx = 0;
  61                                         fprintf( stderr, "jtt %d\n", pamN );
  62                                         --argc;
  63                                         goto nextoption;
  64                                 case 'm':
  65                                         fmodel = 1;
  66                                         break;
  67                                 case 'r':
  68                                         fmodel = -1;
  69                                         break;
  70                                 case 'D':
  71                                         dorp = 'd';
  72                                         break;
  73                                 case 'P':
  74                                         dorp = 'p';
  75                                         break;
  76                 default:
  77                     fprintf( stderr, "illegal option %c\n", c );
  78                     argc = 0;
  79                     break;
  80             }
  81                 }
  82                 nextoption:
  83                         ;
  84         }
  85     if( argc == 1 )
  86     {
  87         cut = atof( (*argv) );
  88         argc--;
  89     }
  90     if( argc != 0 )
  91     {
  92         fprintf( stderr, "options: Check source file !\n" );
  93         exit( 1 );
  94     }
  95         if( tbitr == 1 && outgap == 0 )
  96         {
  97                 fprintf( stderr, "conflicting options : o, m or u\n" );
  98                 exit( 1 );
  99         }
 100         if( alg == 'C' && outgap == 0 )
 101         {
 102                 fprintf( stderr, "conflicting options : C, o\n" );
 103                 exit( 1 );
 104         }
 105 }
 106
 107 int main( int argc, char *argv[] )
 108 {
 109         static int  nlen[M];
 110         static char name[M][B], **seq;
 111         int i, j, alloclen, c;
 112         double **mtx;
 113         double *self;
 114         double tmpdouble;
 115         FILE *fp;
 116
 117         arguments( argc, argv );
 118
 119         getnumlen( stdin );
 120         rewind( stdin );
 121
 122         if( njob < 2 )
 123         {
 124                 fprintf( stderr, "At least 2 sequences should be input!\n"
 125                                                  "Only %d sequence found.\n", njob );
 126                 exit( 1 );
 127         }
 128
 129         seq = AllocateCharMtx( njob, nlenmax*9+1 );
 130         mtx = AllocateDoubleMtx( njob, njob );
 131         self = AllocateDoubleVec( njob );
 132         alloclen = nlenmax*9;
 133
 134         readData( stdin, name, nlen, seq );
 135         constants( njob, seq );
 136
 137
 138
 139
 140         c = seqcheck( seq );
 141         if( c )
 142         {
 143                 fprintf( stderr, "Illeagal character %c\n", c );
 144                 exit( 1 );
 145         }
 146
 147         for( i=0; i<njob; i++ )
 148         {
 149                 self[i] = (double)substitution_nid( seq[i], seq[i] );
 150 //              fprintf( stdout, "self[%d] = %f\n", i, self[i] );
 151         }
 152
 153         for( i=0; i<njob-1; i++ )
 154                 for( j=i+1; j<njob; j++ )
 155                 {
 156                         tmpdouble = (double)substitution_score( seq[i], seq[j] );
 157 //                      fprintf( stdout, "tmpdouble = %f\n", tmpdouble );
 158                         mtx[i][j] = ( 1.0 - tmpdouble / MIN( self[i], self[j] ) );
 159                         if( mtx[i][j] < 0.95 )
 160                                 mtx[i][j] = - log( 1.0 - mtx[i][j] );
 161                         else
 162                                 mtx[i][j] = 3.0;
 163                 }
 164
 165 #if TEST
 166         for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 167                 fprintf( stdout, "i=%d, j=%d, mtx[][] = %f\n", i, j, mtx[i][j] );
 168 #endif
 169
 170         fp = fopen( "hat2", "w" );
 171         WriteHat2( fp, njob, name, mtx );
 172         fclose( fp );
 173         exit( 0 );
 174
 175         return( 0 );
 176 }