source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/mccaskillwrap.c

Last change on this file was 10371, checked in by aboeckma, 11 years ago

updated mafft version. Added extensions (no svn ignore, yet)
File size: 10.5 KB
Line 
1#include "mltaln.h"
2
3#define DEBUG 0
4
5static char *whereismccaskillmea;
6
7#ifdef enablemultithread
8typedef 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
21void 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
33static 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
77static 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                        break;
106                }
107                *jobpospt = i+1;
108                pthread_mutex_unlock( targ->mutex );
109
110                commongappick_record( 1, nogap+i, gapmap[i] );
111                if( strlen( nogap[i] ) == 0 ) continue;
112
113                sprintf( dirname, "_%d", i );
114                sprintf( com, "rm -rf %s", dirname );
115                system( com );
116                sprintf( com, "mkdir %s", dirname );
117                system( com );
118
119                fprintf( stderr, "%d / %d (by thread %4d)\n", i+1, njob, thread_no );
120                sprintf( com, "%s/_mccaskillinorg", dirname );
121                infp = fopen( com, "w" );
122//              fprintf( infp, ">in\n%s\n", nogap[i] );
123                fprintf( infp, ">in\n" );
124                write1seq( infp, nogap[i] );
125                fclose( infp );
126
127                sprintf( com, "tr -d '\\r' < %s/_mccaskillinorg > %s/_mccaskillin", dirname, dirname );
128                system( com ); // for cygwin, wakaran
129                if( alg == 'G' )
130                        sprintf( com, "cd %s; %s/dafs --mafft-out _mccaskillout _mccaskillin > _dum1 2>_dum", dirname, whereismccaskillmea );
131                else
132                        sprintf( com, "cd %s; %s/mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", dirname, whereismccaskillmea );
133                res = system( com );
134
135                if( res )
136                {
137                        fprintf( stderr, "ERROR IN mccaskill_mea\n" );
138                        exit( 1 );
139                }
140
141                sprintf( com, "%s/_mccaskillout", dirname );
142                infp = fopen( com, "r" );
143                readrawmccaskill( infp, pairprob[i], nlenmax );
144                fclose( infp );
145
146                sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
147                if( system( com ) )
148                {
149                        fprintf( stderr, "retrying to rmdir\n" );
150//                      nanosleep( 100000 );
151                        sleep( 1 );
152                        system( com );
153                }
154        }
155        free( dirname );
156        free( com );
157        return( NULL );
158}
159#endif
160
161void arguments( int argc, char *argv[] )
162{
163    int c;
164        nthread = 1;
165        inputfile = NULL;
166        dorp = NOTSPECIFIED;
167        kimuraR = NOTSPECIFIED;
168        pamN = NOTSPECIFIED;
169        whereismccaskillmea = NULL;
170        alg = 's';
171
172    while( --argc > 0 && (*++argv)[0] == '-' )
173        {
174        while ( (c = *++argv[0]) )
175                {
176            switch( c )
177            {
178                                case 'i':
179                                        inputfile = *++argv;
180                                        fprintf( stderr, "inputfile = %s\n", inputfile );
181                                        --argc;
182                                        goto nextoption;
183                                case 'd':
184                                        whereismccaskillmea = *++argv;
185                                        fprintf( stderr, "whereismccaskillmea = %s\n", whereismccaskillmea );
186                                        --argc;
187                                        goto nextoption;
188                                case 'C':
189                                        nthread = atoi( *++argv );
190                                        fprintf( stderr, "nthread = %d\n", nthread );
191                                        --argc; 
192                                        goto nextoption;
193                                case 's':
194                                        alg = 's'; // use scarna; default
195                                        break;
196                                case 'G':
197                                        alg = 'G'; // use dafs, instead of scarna
198                                        break;
199                default:
200                    fprintf( stderr, "illegal option %c\n", c );
201                    argc = 0;
202                    break;
203            }
204                }
205                nextoption:
206                        ;
207        }
208    if( argc != 0 ) 
209    {
210        fprintf( stderr, "options: Check source file !\n" );
211        exit( 1 );
212    }
213}
214
215
216int main( int argc, char *argv[] )
217{
218        static char com[10000];
219        static int  *nlen;     
220        int left, right;
221        int res;
222        static char **name, **seq, **nogap;
223        static int **gapmap;
224        static int *order;
225        int i, j;
226        FILE *infp;
227        RNApair ***pairprob;
228        RNApair **alnpairprob;
229        RNApair *pairprobpt;
230        RNApair *pt;
231        int *alnpairnum;
232        float prob;
233        int adpos;
234
235        arguments( argc, argv );
236#ifndef enablemultithread
237        nthread = 0;
238#endif
239
240        if( inputfile )
241        {
242                infp = fopen( inputfile, "r" );
243                if( !infp )
244                {
245                        fprintf( stderr, "Cannot open %s\n", inputfile );
246                        exit( 1 );
247                }
248        }
249        else
250                infp = stdin;
251
252        if( !whereismccaskillmea )
253                whereismccaskillmea = "";
254
255        getnumlen( infp );
256        rewind( infp );
257
258        if( dorp != 'd' )
259        {
260                fprintf( stderr, "nuc only\n" );
261                exit( 1 );
262        }
263
264        seq = AllocateCharMtx( njob, nlenmax*2+1 );
265        nogap = AllocateCharMtx( njob, nlenmax*2+1 );
266        gapmap = AllocateIntMtx( njob, nlenmax*2+1 );
267        order = AllocateIntVec( njob );
268        name = AllocateCharMtx( njob, B+1 );
269        nlen = AllocateIntVec( njob );
270        pairprob = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
271        alnpairprob = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
272        alnpairnum = AllocateIntVec( nlenmax );
273
274        for( i=0; i<nlenmax; i++ ) alnpairnum[i] = 0;
275
276        readData_pointer( infp, name, nlen, seq );
277        fclose( infp );
278
279        for( i=0; i<njob; i++ )
280        {
281                pairprob[i] = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
282                for( j=0; j<nlenmax; j++ )
283                {
284                        pairprob[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
285                        pairprob[i][j][0].bestpos = -1;
286                        pairprob[i][j][0].bestscore = -1.0;
287                }
288                strcpy( nogap[i], seq[i] );
289                order[i] = i;
290        }
291        for( j=0; j<nlenmax; j++ )
292        {
293                alnpairprob[j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
294                alnpairprob[j][0].bestpos = -1;
295                alnpairprob[j][0].bestscore = -1.0;
296        }
297
298
299        constants( njob, seq );
300
301        if( alg == 'G' )
302                fprintf( stderr, "Running DAFS (Sato et al. 2012; http://www.ncrna.org/).\n" );
303        else
304                fprintf( stderr, "Running mxscarna with the mccaskill_mea mode.\n" );
305#ifdef enablemultithread
306        if( nthread > 0 )
307        {
308                int jobpos;
309                pthread_t *handle;
310                pthread_mutex_t mutex;
311                thread_arg_t *targ;
312                jobpos = 0;
313
314                targ = calloc( nthread, sizeof( thread_arg_t ) );
315                handle = calloc( nthread, sizeof( pthread_t ) );
316                pthread_mutex_init( &mutex, NULL );
317
318                for( i=0; i<nthread; i++ )
319                {
320                        targ[i].thread_no = i;
321                        targ[i].njob = njob;
322                        targ[i].jobpospt = &jobpos;
323                        targ[i].gapmap = gapmap;
324                        targ[i].nogap = nogap;
325                        targ[i].nlenmax = nlenmax;
326                        targ[i].pairprob = pairprob;
327                        targ[i].mutex = &mutex;
328
329//                      athread( targ );
330                        pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
331                       
332                }
333
334                for( i=0; i<nthread; i++ )
335                {
336                        pthread_join( handle[i], NULL );
337                }
338                pthread_mutex_destroy( &mutex );
339
340                free( handle );
341                free( targ );
342
343
344                for( i=0; i<njob; i++ )
345                {
346                        fprintf( stdout, ">%d\n", i );
347                        outmccaskill( stdout, pairprob[i], nlenmax );
348                }
349        }
350        else
351#endif
352        {
353                for( i=0; i<njob; i++ )
354                {
355                        fprintf( stderr, "%d / %d\n", i+1, njob );
356                        commongappick_record( 1, nogap+i, gapmap[i] );
357                        if( strlen( nogap[i] ) == 0 ) 
358                        {
359                                fprintf( stdout, ">%d\n", i );
360                                continue;
361                        }
362
363                        infp = fopen( "_mccaskillinorg", "w" );
364//                      fprintf( infp, ">in\n%s\n", nogap[i] );
365                        fprintf( infp, ">in\n" );
366                        write1seq( infp, nogap[i] );
367                        fclose( infp );
368       
369                        system( "tr -d '\\r' < _mccaskillinorg > _mccaskillin" ); // for cygwin, wakaran
370                        if( alg == 'G' )
371                                sprintf( com, "env PATH=%s dafs --mafft-out _mccaskillout _mccaskillin > _dum1 2>_dum", whereismccaskillmea );
372                        else
373                                sprintf( com, "env PATH=%s mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", whereismccaskillmea );
374                        res = system( com );
375       
376                        if( res )
377                        {
378                                fprintf( stderr, "ERROR IN mccaskill_mea\n" );
379                                exit( 1 );
380                        }
381       
382                        infp = fopen( "_mccaskillout", "r" );
383                        readrawmccaskill( infp, pairprob[i], nlenmax );
384                        fclose( infp );
385                        fprintf( stdout, ">%d\n", i );
386                        outmccaskill( stdout, pairprob[i], nlenmax );
387                }
388        }
389
390        for( i=0; i<njob; i++ )
391        {
392                for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
393                {
394                        left = gapmap[i][j];
395                        right = gapmap[i][pairprobpt->bestpos];
396                        prob = pairprobpt->bestscore;
397
398                        for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
399                                if( pt->bestpos == right ) break;
400
401                        if( pt->bestpos == -1 )
402                        {
403                                alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
404                                adpos = alnpairnum[left];
405                                alnpairnum[left]++;
406                                alnpairprob[left][adpos].bestscore = 0.0;
407                                alnpairprob[left][adpos].bestpos = right;
408                                alnpairprob[left][adpos+1].bestscore = -1.0;
409                                alnpairprob[left][adpos+1].bestpos = -1;
410                                pt = alnpairprob[left]+adpos;
411                        }
412                        else
413                                adpos = pt-alnpairprob[left];
414
415                        pt->bestscore += prob;
416                        if( pt->bestpos != right )
417                        {
418                                fprintf( stderr, "okashii!\n" );
419                                exit( 1 );
420                        }
421//                      fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
422                }
423        }
424
425        for( i=0; i<njob; i++ )
426        {
427                for( j=0; j<nlenmax; j++ ) free( pairprob[i][j] );
428                free( pairprob[i] );
429        }
430        free( pairprob );
431        for( j=0; j<nlenmax; j++ ) free( alnpairprob[j] );
432        free( alnpairprob );
433        free( alnpairnum );
434        fprintf( stderr, "%d thread(s)\n", nthread );
435        return( 0 );
436
437#if 0
438        fprintf( stdout, "result=\n" );
439
440        for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
441        {
442                pairprobpt->bestscore /= (float)njob;
443                left = i;
444                right = pairprobpt->bestpos;
445                prob = pairprobpt->bestscore;
446                fprintf( stdout, "%d-%d, %f\n", left, right, prob );
447        }
448
449        return( 0 );
450#endif
451}
Note: See TracBrowser for help on using the repository browser.