5 static char *whereismccaskillmea;
7 #ifdef enablemultithread
8 typedef struct _thread_arg
17 pthread_mutex_t *mutex;
21 void outmccaskill( FILE *fp, RNApair **pairprob, int length )
25 for( i=0; i<length; i++ ) for( pt=pairprob[i]; pt->bestpos!=-1; pt++ )
28 fprintf( fp, "%d %d %50.40f\n", i, pt->bestpos, pt->bestscore );
33 static void readrawmccaskill( FILE *fp, RNApair **pairprob, int length )
41 pairnum = (int *)calloc( length, sizeof( int ) );
42 for( i=0; i<length; i++ ) pairnum[i] = 0;
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 );
53 if( left != right && prob > 0.0 )
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;
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 );
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;
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 );
76 #ifdef enablemultithread
77 static void *athread( void *arg )
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;
93 dirname = calloc( 100, sizeof( char ) );
94 com = calloc( 1000, sizeof( char ) );
99 pthread_mutex_lock( targ->mutex );
103 pthread_mutex_unlock( targ->mutex );
107 pthread_mutex_unlock( targ->mutex );
110 sprintf( dirname, "_%d", i );
111 sprintf( com, "rm -rf %s", dirname );
113 sprintf( com, "mkdir %s", dirname );
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] );
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 );
132 fprintf( stderr, "ERROR IN mccaskill_mea\n" );
136 sprintf( com, "%s/_mccaskillout", dirname );
137 infp = fopen( com, "r" );
138 readrawmccaskill( infp, pairprob[i], nlenmax );
141 sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
144 fprintf( stderr, "retrying to rmdir\n" );
145 // nanosleep( 100000 );
155 void arguments( int argc, char *argv[] )
161 kimuraR = NOTSPECIFIED;
163 whereismccaskillmea = NULL;
165 while( --argc > 0 && (*++argv)[0] == '-' )
167 while ( (c = *++argv[0]) )
173 fprintf( stderr, "inputfile = %s\n", inputfile );
177 whereismccaskillmea = *++argv;
178 fprintf( stderr, "whereismccaskillmea = %s\n", whereismccaskillmea );
182 nthread = atoi( *++argv );
183 fprintf( stderr, "nthread = %d\n", nthread );
187 fprintf( stderr, "illegal option %c\n", c );
197 fprintf( stderr, "options: Check source file !\n" );
203 int main( int argc, char *argv[] )
205 static char com[10000];
209 static char **name, **seq, **nogap;
215 RNApair **alnpairprob;
222 arguments( argc, argv );
223 #ifndef enablemultithread
229 infp = fopen( inputfile, "r" );
232 fprintf( stderr, "Cannot open %s\n", inputfile );
239 if( !whereismccaskillmea )
240 whereismccaskillmea = "";
247 fprintf( stderr, "nuc only\n" );
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 );
261 for( i=0; i<nlenmax; i++ ) alnpairnum[i] = 0;
263 readData_pointer( infp, name, nlen, seq );
266 for( i=0; i<njob; i++ )
268 pairprob[i] = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
269 for( j=0; j<nlenmax; j++ )
271 pairprob[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
272 pairprob[i][j][0].bestpos = -1;
273 pairprob[i][j][0].bestscore = -1.0;
275 strcpy( nogap[i], seq[i] );
278 for( j=0; j<nlenmax; j++ )
280 alnpairprob[j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
281 alnpairprob[j][0].bestpos = -1;
282 alnpairprob[j][0].bestscore = -1.0;
286 constants( njob, seq );
288 fprintf( stderr, "Running mxscarna with the mccaskill_mea mode.\n" );
289 #ifdef enablemultithread
294 pthread_mutex_t mutex;
298 targ = calloc( nthread, sizeof( thread_arg_t ) );
299 handle = calloc( nthread, sizeof( pthread_t ) );
300 pthread_mutex_init( &mutex, NULL );
302 for( i=0; i<nthread; i++ )
304 targ[i].thread_no = i;
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;
314 pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
318 for( i=0; i<nthread; i++ )
320 pthread_join( handle[i], NULL );
322 pthread_mutex_destroy( &mutex );
325 for( i=0; i<njob; i++ )
327 fprintf( stdout, ">%d\n", i );
328 outmccaskill( stdout, pairprob[i], nlenmax );
334 for( i=0; i<njob; i++ )
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] );
344 system( "tr -d '\\r' < _mccaskillinorg > _mccaskillin" ); // for cygwin, wakaran
345 sprintf( com, "env PATH=%s mxscarnamod -m -writebpp _mccaskillin > _mccaskillout 2>_dum", whereismccaskillmea );
350 fprintf( stderr, "ERROR IN mccaskill_mea\n" );
354 infp = fopen( "_mccaskillout", "r" );
355 readrawmccaskill( infp, pairprob[i], nlenmax );
357 fprintf( stdout, ">%d\n", i );
358 outmccaskill( stdout, pairprob[i], nlenmax );
362 for( i=0; i<njob; i++ )
364 for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
367 right = gapmap[i][pairprobpt->bestpos];
368 prob = pairprobpt->bestscore;
370 for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
371 if( pt->bestpos == right ) break;
373 if( pt->bestpos == -1 )
375 alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
376 adpos = 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;
385 adpos = pt-alnpairprob[left];
387 pt->bestscore += prob;
388 if( pt->bestpos != right )
390 fprintf( stderr, "okashii!\n" );
393 // fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
397 for( i=0; i<njob; i++ )
399 for( j=0; j<nlenmax; j++ ) free( pairprob[i][j] );
403 for( j=0; j<nlenmax; j++ ) free( alnpairprob[j] );
406 fprintf( stderr, "%d thread(s)\n", nthread );
410 fprintf( stdout, "result=\n" );
412 for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
414 pairprobpt->bestscore /= (float)njob;
416 right = pairprobpt->bestpos;
417 prob = pairprobpt->bestscore;
418 fprintf( stdout, "%d-%d, %f\n", left, right, prob );