binaries/src/mafft/core/multi2hat3s.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 1
   6 #define TSUYOSAFACTOR 100
   7
   8
   9 static int nhomologs;
  10 static int seedoffset;
  11
  12 void strip( char *s )
  13 {
  14         char *pt = s;
  15         while( *++pt )
  16                 if( *pt == '\n' ) *pt = 0;
  17 }
  18
  19
  20 void arguments( int argc, char *argv[] )
  21 {
  22     int c;
  23
  24         seedoffset = 0;
  25         nhomologs = 1;
  26         inputfile = NULL;
  27         fftkeika = 0;
  28         pslocal = -1000.0;
  29         constraint = 0;
  30         nblosum = 62;
  31         fmodel = 0;
  32         calledByXced = 0;
  33         devide = 0;
  34         use_fft = 0;
  35         fftscore = 1;
  36         fftRepeatStop = 0;
  37         fftNoAnchStop = 0;
  38     weight = 3;
  39     utree = 1;
  40         tbutree = 1;
  41     refine = 0;
  42     check = 1;
  43     cut = 0.0;
  44     disp = 0;
  45     outgap = 1;
  46     alg = 'A';
  47     mix = 0;
  48         tbitr = 0;
  49         scmtd = 5;
  50         tbweight = 0;
  51         tbrweight = 3;
  52         checkC = 0;
  53         treemethod = 'x';
  54         contin = 0;
  55         scoremtx = 1;
  56         kobetsubunkatsu = 0;
  57         divpairscore = 0;
  58         dorp = NOTSPECIFIED;
  59         ppenalty = NOTSPECIFIED;
  60         ppenalty_OP = NOTSPECIFIED;
  61         ppenalty_ex = NOTSPECIFIED;
  62         ppenalty_EX = NOTSPECIFIED;
  63         poffset = NOTSPECIFIED;
  64         kimuraR = NOTSPECIFIED;
  65         pamN = NOTSPECIFIED;
  66         geta2 = GETA2;
  67         fftWinSize = NOTSPECIFIED;
  68         fftThreshold = NOTSPECIFIED;
  69
  70     while( --argc > 0 && (*++argv)[0] == '-' )
  71         {
  72         while ( ( c = *++argv[0] ) )
  73                 {
  74             switch( c )
  75             {
  76                                 case 'i':
  77                                         inputfile = *++argv;
  78                                         fprintf( stderr, "seed = %s\n", inputfile );
  79                                         --argc;
  80                                         goto nextoption;
  81                                 case 't':
  82                                         nhomologs = atoi( *++argv );
  83                                         fprintf( stderr, "nhomologs = %d\n", nhomologs );
  84                                         --argc;
  85                                         goto nextoption;
  86                                 case 'o':
  87                                         seedoffset = atoi( *++argv );
  88                                         fprintf( stderr, "seedoffset = %d\n", seedoffset );
  89                                         --argc;
  90                                         goto nextoption;
  91                                 case 'D':
  92                                         dorp = 'd';
  93                                         break;
  94                                 case 'P':
  95                                         dorp = 'p';
  96                                         break;
  97                 default:
  98                     fprintf( stderr, "illegal option %c\n", c );
  99                     argc = 0;
 100                     break;
 101             }
 102                 }
 103                 nextoption:
 104                         ;
 105         }
 106     if( argc == 1 )
 107     {
 108         cut = atof( (*argv) );
 109         argc--;
 110     }
 111     if( argc != 0 )
 112     {
 113         fprintf( stderr, "options: Check source file !\n" );
 114         exit( 1 );
 115     }
 116         if( tbitr == 1 && outgap == 0 )
 117         {
 118                 fprintf( stderr, "conflicting options : o, m or u\n" );
 119                 exit( 1 );
 120         }
 121         if( alg == 'C' && outgap == 0 )
 122         {
 123                 fprintf( stderr, "conflicting options : C, o\n" );
 124                 exit( 1 );
 125         }
 126 }
 127
 128 int countamino( char *s, int end )
 129 {
 130         int val = 0;
 131         while( end-- )
 132                 if( *s++ != '-' ) val++;
 133         return( val );
 134 }
 135
 136 static void pairalign( char name[M][B], int nlen[M], char **seq, double *effarr, int alloclen )
 137 {
 138         int i, j;
 139         FILE *hat3p;
 140         float pscore = 0.0; // by D.Mathog
 141         static double *effarr1 = NULL;
 142         static double *effarr2 = NULL;
 143         char *aseq;
 144         static char **pseq;
 145         LocalHom **localhomtable, *tmpptr;
 146         double tsuyosa;
 147
 148         if( nhomologs < 1 ) nhomologs = 1; // tsuyosa=0.0 wo sakeru
 149         tsuyosa = (double)nhomologs * nhomologs * TSUYOSAFACTOR;
 150         fprintf( stderr, "tsuyosa = %f\n", tsuyosa );
 151         localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 152         for( i=0; i<njob; i++)
 153         {
 154                 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
 155                 for( j=0; j<njob; j++)
 156                 {
 157                         localhomtable[i][j].start1 = -1;
 158                         localhomtable[i][j].end1 = -1;
 159                         localhomtable[i][j].start2 = -1;
 160                         localhomtable[i][j].end2 = -1;
 161                         localhomtable[i][j].opt = -1.0;
 162                         localhomtable[i][j].next = NULL;
 163                 }
 164         }
 165
 166         if( effarr1 == NULL )
 167         {
 168                 effarr1 = AllocateDoubleVec( njob );
 169                 effarr2 = AllocateDoubleVec( njob );
 170                 pseq = AllocateCharMtx( 2, 0 );
 171                 aseq = AllocateCharVec( nlenmax*9+1 );
 172 #if 0
 173 #else
 174 #endif
 175         }
 176
 177 #if 0
 178         fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
 179 #endif
 180
 181 #if 0
 182         for( i=0; i<njob; i++ )
 183                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
 184 #endif
 185
 186
 187 //      writePre( njob, name, nlen, aseq, 0 );
 188
 189         hat3p = fopen( "hat3", "w" );
 190         if( !hat3p ) ErrorExit( "Cannot open hat3." );
 191         for( i=0; i<njob-1; i++ )
 192         {
 193                 for( j=i+1; j<njob; j++ )
 194                 {
 195                         pseq[0] = seq[i];
 196                         pseq[1] = seq[j];
 197
 198                         if( strlen( pseq[0] ) != strlen( pseq[1] ) )
 199                         {
 200                                 fprintf( stderr, "## ERROR  ###\n" );
 201                                 fprintf( stderr, "Not aligned,  %s - %s\n", name[i], name[j] );
 202                                 fprintf( stderr, "## ERROR  ###\n" );
 203                                 exit( 1 );
 204                         }
 205
 206
 207                         fprintf( stderr, "adding %d-%d\n", i, j );
 208                         putlocalhom2( pseq[0], pseq[1], localhomtable[i]+j, 0, 0, (int)pscore, strlen( pseq[0] ) );
 209                         for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
 210                         {
 211                                 if( tmpptr->opt == -1.0 ) continue;
 212                                 if( tmpptr->start1 == -1 ) continue;
 213                                 fprintf( hat3p, "%d %d %d %6.3f %d %d %d %d k\n", i+seedoffset, j+seedoffset, tmpptr->overlapaa, tmpptr->opt * tsuyosa, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2 );
 214                         }
 215                 }
 216         }
 217         fclose( hat3p );
 218
 219 #if DEBUG
 220         fprintf( stderr, "calling FreeLocalHomTable\n" );
 221 #endif
 222         FreeLocalHomTable( localhomtable, njob );
 223 #if DEBUG
 224         fprintf( stderr, "done. FreeLocalHomTable\n" );
 225 #endif
 226 }
 227
 228 static void WriteOptions( FILE *fp )
 229 {
 230
 231         if( dorp == 'd' ) fprintf( fp, "DNA\n" );
 232         else
 233         {
 234                 if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
 235                 else if( scoremtx ==  1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
 236                 else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
 237         }
 238     fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 239     if( use_fft ) fprintf( fp, "FFT on\n" );
 240
 241         fprintf( fp, "tree-base method\n" );
 242         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
 243         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
 244         if( tbitr || tbweight )
 245         {
 246                 fprintf( fp, "iterate at each step\n" );
 247                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
 248                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
 249                 if( tbweight ) fprintf( fp, "  weighted\n" );
 250                 fprintf( fp, "\n" );
 251         }
 252
 253          fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
 254
 255         if( alg == 'a' )
 256                 fprintf( fp, "Algorithm A\n" );
 257         else if( alg == 'A' )
 258                 fprintf( fp, "Algorithm A+\n" );
 259         else if( alg == 'S' )
 260                 fprintf( fp, "Apgorithm S\n" );
 261         else if( alg == 'C' )
 262                 fprintf( fp, "Apgorithm A+/C\n" );
 263         else
 264                 fprintf( fp, "Unknown algorithm\n" );
 265
 266         if( treemethod == 'x' )
 267                 fprintf( fp, "Tree = UPGMA (3).\n" );
 268         else if( treemethod == 's' )
 269                 fprintf( fp, "Tree = UPGMA (2).\n" );
 270         else if( treemethod == 'p' )
 271                 fprintf( fp, "Tree = UPGMA (1).\n" );
 272         else
 273                 fprintf( fp, "Unknown tree.\n" );
 274
 275     if( use_fft )
 276     {
 277         fprintf( fp, "FFT on\n" );
 278         if( dorp == 'd' )
 279             fprintf( fp, "Basis : 4 nucleotides\n" );
 280         else
 281         {
 282             if( fftscore )
 283                 fprintf( fp, "Basis : Polarity and Volume\n" );
 284             else
 285                 fprintf( fp, "Basis : 20 amino acids\n" );
 286         }
 287         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
 288         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
 289     }
 290         else
 291         fprintf( fp, "FFT off\n" );
 292         fflush( fp );
 293 }
 294
 295
 296 int main( int argc, char *argv[] )
 297 {
 298         static int  nlen[M];
 299         static char name[M][B], **seq;
 300         static char **bseq;
 301         static double *eff;
 302         int i;
 303         char c;
 304         int alloclen;
 305         FILE *infp;
 306
 307         arguments( argc, argv );
 308
 309         if( inputfile )
 310         {
 311                 infp = fopen( inputfile, "r" );
 312                 if( !infp )
 313                 {
 314                         fprintf( stderr, "Cannot open %s\n", inputfile );
 315                         exit( 1 );
 316                 }
 317         }
 318         else
 319                 infp = stdin;
 320
 321         getnumlen( infp );
 322         rewind( infp );
 323
 324         if( njob < 2 )
 325         {
 326                 fprintf( stderr, "At least 2 sequences should be input!\n"
 327                                                  "Only %d sequence found.\n", njob );
 328                 exit( 1 );
 329         }
 330
 331         seq = AllocateCharMtx( njob, nlenmax*9+1 );
 332         bseq = AllocateCharMtx( njob, nlenmax*9+1 );
 333         alloclen = nlenmax*9;
 334
 335         eff = AllocateDoubleVec( njob );
 336
 337 #if 0
 338         Read( name, nlen, seq );
 339 #else
 340         readData( infp, name, nlen, seq );
 341 #endif
 342         fclose( infp );
 343
 344         constants( njob, seq );
 345
 346 #if 0
 347         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
 348 #endif
 349
 350         initSignalSM();
 351
 352         initFiles();
 353
 354         WriteOptions( trap_g );
 355
 356         c = seqcheck( seq );
 357         if( c )
 358         {
 359                 fprintf( stderr, "Illeagal character %c\n", c );
 360                 exit( 1 );
 361         }
 362
 363 //      writePre( njob, name, nlen, seq, 0 );
 364
 365         for( i=0; i<njob; i++ ) eff[i] = 1.0;
 366
 367
 368         for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
 369
 370
 371         for( i=0; i<njob; i++ ) fprintf( stdout, ">_seed_%s\n%s\n", name[i]+1, bseq[i] );
 372
 373         pairalign( name, nlen, seq, eff, alloclen );
 374
 375         fprintf( trap_g, "done.\n" );
 376 #if DEBUG
 377         fprintf( stderr, "closing trap_g\n" );
 378 #endif
 379         fclose( trap_g );
 380
 381 #if IODEBUG
 382         fprintf( stderr, "OSHIMAI\n" );
 383 #endif
 384         SHOWVERSION;
 385         return( 0 );
 386 }