website/archive/binaries/mac/src/mafft/core/mccaskillwrap.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4
   5 static char *whereismccaskillmea;
   6
   7 #ifdef enablemultithread
   8 typedef struct _thread_arg
   9 {
  10         int thread_no;
  11         int njob;
  12         int *jobpospt;
  13         int **gapmap;
  14         char **nogap;
  15         int nlenmax;
  16         RNApair ***pairprob;
  17         pthread_mutex_t *mutex;
  18 } thread_arg_t;
  19 #endif
  20
  21 void outmccaskill( FILE *fp, RNApair **pairprob, int length )
  22 {
  23         int i;
  24         RNApair *pt;
  25         for( i=0; i<length; i++ ) for( pt=pairprob[i]; pt->bestpos!=-1; pt++ )
  26         {
  27                 if( pt->bestpos > i )
  28                         fprintf( fp, "%d %d %50.40f\n", i, pt->bestpos, pt->bestscore );
  29         }
  30 }
  31
  32 #if 1
  33 static void readrawmccaskill( FILE *fp, RNApair **pairprob, int length )
  34 {
  35         char gett[1000];
  36         int *pairnum;
  37         int i;
  38         int left, right;
  39         float prob;
  40
  41         pairnum = (int *)calloc( length, sizeof( int ) );
  42         for( i=0; i<length; i++ ) pairnum[i] = 0;
  43
  44         while( 1 )
  45         {
  46                 fgets( gett, 999, fp );
  47                 if( feof( fp ) ) break;
  48                 if( gett[0] == '>' ) continue;
  49                 sscanf( gett, "%d %d %f", &left, &right, &prob );
  50                 if( prob < 0.01 ) continue; // mxscarna to mafft ryoho ni eikyou
  51 //fprintf( stderr, "gett = %s\n", gett );
  52
  53                 if( left != right && prob > 0.0 )
  54                 {
  55                         pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
  56                         pairprob[left][pairnum[left]].bestscore = prob;
  57                         pairprob[left][pairnum[left]].bestpos = right;
  58                         pairnum[left]++;
  59                         pairprob[left][pairnum[left]].bestscore = -1.0;
  60                         pairprob[left][pairnum[left]].bestpos = -1;
  61 //                      fprintf( stderr, "%d-%d, %f\n", left, right, prob );
  62
  63                         pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
  64                         pairprob[right][pairnum[right]].bestscore = prob;
  65                         pairprob[right][pairnum[right]].bestpos = left;
  66                         pairnum[right]++;
  67                         pairprob[right][pairnum[right]].bestscore = -1.0;
  68                         pairprob[right][pairnum[right]].bestpos = -1;
  69 //                      fprintf( stderr, "%d-%d, %f\n", right, left, prob );
  70                 }
  71         }
  72         free( pairnum );
  73 }
  74 #endif
  75
  76 #ifdef enablemultithread
  77 static void *athread( void *arg )
  78 {
  79         thread_arg_t *targ = (thread_arg_t *)arg;
  80         int thread_no = targ->thread_no;
  81         int njob = targ->njob;
  82         int *jobpospt = targ->jobpospt;
  83         int **gapmap = targ->gapmap;
  84         char **nogap = targ->nogap;
  85         int nlenmax = targ->nlenmax;
  86         RNApair ***pairprob = targ->pairprob;
  87
  88         int i, res;
  89         FILE *infp;
  90         char *com;
  91         char *dirname;
  92
  93         dirname = calloc( 100, sizeof( char ) );
  94         com = calloc( 1000, sizeof( char ) );
  95
  96
  97         while( 1 )
  98         {
  99                 pthread_mutex_lock( targ->mutex );
 100                 i = *jobpospt;
 101                 if( i == njob )
 102                 {
 103                         pthread_mutex_unlock( targ->mutex );
 104                         return( NULL );
 105                 }
 106                 *jobpospt = i+1;
 107                 pthread_mutex_unlock( targ->mutex );
 108
 109
 110                 sprintf( dirname, "_%d", i );
 111                 sprintf( com, "rm -rf %s", dirname );
 112                 system( com );
 113                 sprintf( com, "mkdir %s", dirname );
 114                 system( com );
 115
 116                 fprintf( stderr, "%d / %d (by thread %4d)\n", i+1, njob, thread_no );
 117                 commongappick_record( 1, nogap+i, gapmap[i] );
 118                 sprintf( com, "%s/_mccaskillinorg", dirname );
 119                 infp = fopen( com, "w" );
 120 //              fprintf( infp, ">in\n%s\n", nogap[i] );
 121                 fprintf( infp, ">in\n" );
 122                 write1seq( infp, nogap[i] );
 123                 fclose( infp );
 124
 125                 sprintf( com, "tr -d '\\r' < %s/_mccaskillinorg > %s/_mccaskillin", dirname, dirname );
 126                 system( com ); // for cygwin, wakaran
 127                 sprintf( com, "cd %s; %s/mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", dirname, whereismccaskillmea );
 128                 res = system( com );
 129
 130                 if( res )
 131                 {
 132                         fprintf( stderr, "ERROR IN mccaskill_mea\n" );
 133                         exit( 1 );
 134                 }
 135
 136                 sprintf( com, "%s/_mccaskillout", dirname );
 137                 infp = fopen( com, "r" );
 138                 readrawmccaskill( infp, pairprob[i], nlenmax );
 139                 fclose( infp );
 140
 141                 sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
 142                 if( system( com ) )
 143                 {
 144                         fprintf( stderr, "retrying to rmdir\n" );
 145 //                      nanosleep( 100000 );
 146                         sleep( 1 );
 147                         system( com );
 148                 }
 149         }
 150         free( dirname );
 151         free( com );
 152 }
 153 #endif
 154
 155 void arguments( int argc, char *argv[] )
 156 {
 157     int c;
 158         nthread = 1;
 159         inputfile = NULL;
 160         dorp = NOTSPECIFIED;
 161         kimuraR = NOTSPECIFIED;
 162         pamN = NOTSPECIFIED;
 163         whereismccaskillmea = NULL;
 164
 165     while( --argc > 0 && (*++argv)[0] == '-' )
 166         {
 167         while ( (c = *++argv[0]) )
 168                 {
 169             switch( c )
 170             {
 171                                 case 'i':
 172                                         inputfile = *++argv;
 173                                         fprintf( stderr, "inputfile = %s\n", inputfile );
 174                                         --argc;
 175                                         goto nextoption;
 176                                 case 'd':
 177                                         whereismccaskillmea = *++argv;
 178                                         fprintf( stderr, "whereismccaskillmea = %s\n", whereismccaskillmea );
 179                                         --argc;
 180                                         goto nextoption;
 181                                 case 'C':
 182                                         nthread = atoi( *++argv );
 183                                         fprintf( stderr, "nthread = %d\n", nthread );
 184                                         --argc;
 185                                         goto nextoption;
 186                 default:
 187                     fprintf( stderr, "illegal option %c\n", c );
 188                     argc = 0;
 189                     break;
 190             }
 191                 }
 192                 nextoption:
 193                         ;
 194         }
 195     if( argc != 0 )
 196     {
 197         fprintf( stderr, "options: Check source file !\n" );
 198         exit( 1 );
 199     }
 200 }
 201
 202
 203 int main( int argc, char *argv[] )
 204 {
 205         static char com[10000];
 206         static int  *nlen;
 207         int left, right;
 208         int res;
 209         static char **name, **seq, **nogap;
 210         static int **gapmap;
 211         static int *order;
 212         int i, j;
 213         FILE *infp;
 214         RNApair ***pairprob;
 215         RNApair **alnpairprob;
 216         RNApair *pairprobpt;
 217         RNApair *pt;
 218         int *alnpairnum;
 219         float prob;
 220         int adpos;
 221
 222         arguments( argc, argv );
 223 #ifndef enablemultithread
 224         nthread = 0;
 225 #endif
 226
 227         if( inputfile )
 228         {
 229                 infp = fopen( inputfile, "r" );
 230                 if( !infp )
 231                 {
 232                         fprintf( stderr, "Cannot open %s\n", inputfile );
 233                         exit( 1 );
 234                 }
 235         }
 236         else
 237                 infp = stdin;
 238
 239         if( !whereismccaskillmea )
 240                 whereismccaskillmea = "";
 241
 242         getnumlen( infp );
 243         rewind( infp );
 244
 245         if( dorp != 'd' )
 246         {
 247                 fprintf( stderr, "nuc only\n" );
 248                 exit( 1 );
 249         }
 250
 251         seq = AllocateCharMtx( njob, nlenmax*2+1 );
 252         nogap = AllocateCharMtx( njob, nlenmax*2+1 );
 253         gapmap = AllocateIntMtx( njob, nlenmax*2+1 );
 254         order = AllocateIntVec( njob );
 255         name = AllocateCharMtx( njob, B+1 );
 256         nlen = AllocateIntVec( njob );
 257         pairprob = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 258         alnpairprob = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
 259         alnpairnum = AllocateIntVec( nlenmax );
 260
 261         for( i=0; i<nlenmax; i++ ) alnpairnum[i] = 0;
 262
 263         readData_pointer( infp, name, nlen, seq );
 264         fclose( infp );
 265
 266         for( i=0; i<njob; i++ )
 267         {
 268                 pairprob[i] = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
 269                 for( j=0; j<nlenmax; j++ )
 270                 {
 271                         pairprob[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
 272                         pairprob[i][j][0].bestpos = -1;
 273                         pairprob[i][j][0].bestscore = -1.0;
 274                 }
 275                 strcpy( nogap[i], seq[i] );
 276                 order[i] = i;
 277         }
 278         for( j=0; j<nlenmax; j++ )
 279         {
 280                 alnpairprob[j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
 281                 alnpairprob[j][0].bestpos = -1;
 282                 alnpairprob[j][0].bestscore = -1.0;
 283         }
 284
 285
 286         constants( njob, seq );
 287
 288         fprintf( stderr, "Running mxscarna with the mccaskill_mea mode.\n" );
 289 #ifdef enablemultithread
 290         if( nthread > 0 )
 291         {
 292                 int jobpos;
 293                 pthread_t *handle;
 294                 pthread_mutex_t mutex;
 295                 thread_arg_t *targ;
 296                 jobpos = 0;
 297
 298                 targ = calloc( nthread, sizeof( thread_arg_t ) );
 299                 handle = calloc( nthread, sizeof( pthread_t ) );
 300                 pthread_mutex_init( &mutex, NULL );
 301
 302                 for( i=0; i<nthread; i++ )
 303                 {
 304                         targ[i].thread_no = i;
 305                         targ[i].njob = njob;
 306                         targ[i].jobpospt = &jobpos;
 307                         targ[i].gapmap = gapmap;
 308                         targ[i].nogap = nogap;
 309                         targ[i].nlenmax = nlenmax;
 310                         targ[i].pairprob = pairprob;
 311                         targ[i].mutex = &mutex;
 312
 313 //                      athread( targ );
 314                         pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
 315
 316                 }
 317
 318                 for( i=0; i<nthread; i++ )
 319                 {
 320                         pthread_join( handle[i], NULL );
 321                 }
 322                 pthread_mutex_destroy( &mutex );
 323
 324
 325                 for( i=0; i<njob; i++ )
 326                 {
 327                         fprintf( stdout, ">%d\n", i );
 328                         outmccaskill( stdout, pairprob[i], nlenmax );
 329                 }
 330         }
 331         else
 332 #endif
 333         {
 334                 for( i=0; i<njob; i++ )
 335                 {
 336                         fprintf( stderr, "%d / %d\n", i+1, njob );
 337                         commongappick_record( 1, nogap+i, gapmap[i] );
 338                         infp = fopen( "_mccaskillinorg", "w" );
 339 //                      fprintf( infp, ">in\n%s\n", nogap[i] );
 340                         fprintf( infp, ">in\n" );
 341                         write1seq( infp, nogap[i] );
 342                         fclose( infp );
 343
 344                         system( "tr -d '\\r' < _mccaskillinorg > _mccaskillin" ); // for cygwin, wakaran
 345                         sprintf( com, "env PATH=%s mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", whereismccaskillmea );
 346                         res = system( com );
 347
 348                         if( res )
 349                         {
 350                                 fprintf( stderr, "ERROR IN mccaskill_mea\n" );
 351                                 exit( 1 );
 352                         }
 353
 354                         infp = fopen( "_mccaskillout", "r" );
 355                         readrawmccaskill( infp, pairprob[i], nlenmax );
 356                         fclose( infp );
 357                         fprintf( stdout, ">%d\n", i );
 358                         outmccaskill( stdout, pairprob[i], nlenmax );
 359                 }
 360         }
 361
 362         for( i=0; i<njob; i++ )
 363         {
 364                 for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
 365                 {
 366                         left = gapmap[i][j];
 367                         right = gapmap[i][pairprobpt->bestpos];
 368                         prob = pairprobpt->bestscore;
 369
 370                         for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
 371                                 if( pt->bestpos == right ) break;
 372
 373                         if( pt->bestpos == -1 )
 374                         {
 375                                 alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
 376                                 adpos = alnpairnum[left];
 377                                 alnpairnum[left]++;
 378                                 alnpairprob[left][adpos].bestscore = 0.0;
 379                                 alnpairprob[left][adpos].bestpos = right;
 380                                 alnpairprob[left][adpos+1].bestscore = -1.0;
 381                                 alnpairprob[left][adpos+1].bestpos = -1;
 382                                 pt = alnpairprob[left]+adpos;
 383                         }
 384                         else
 385                                 adpos = pt-alnpairprob[left];
 386
 387                         pt->bestscore += prob;
 388                         if( pt->bestpos != right )
 389                         {
 390                                 fprintf( stderr, "okashii!\n" );
 391                                 exit( 1 );
 392                         }
 393 //                      fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
 394                 }
 395         }
 396
 397         for( i=0; i<njob; i++ )
 398         {
 399                 for( j=0; j<nlenmax; j++ ) free( pairprob[i][j] );
 400                 free( pairprob[i] );
 401         }
 402         free( pairprob );
 403         for( j=0; j<nlenmax; j++ ) free( alnpairprob[j] );
 404         free( alnpairprob );
 405         free( alnpairnum );
 406         fprintf( stderr, "%d thread(s)\n", nthread );
 407         return( 0 );
 408
 409 #if 0
 410         fprintf( stdout, "result=\n" );
 411
 412         for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
 413         {
 414                 pairprobpt->bestscore /= (float)njob;
 415                 left = i;
 416                 right = pairprobpt->bestpos;
 417                 prob = pairprobpt->bestscore;
 418                 fprintf( stdout, "%d-%d, %f\n", left, right, prob );
 419         }
 420
 421         return( 0 );
 422 #endif
 423 }