source: branches/items/GDE/MAFFT/mafft-7.055-with-extensions/core/makedirectionlist.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: 17.4 KB
Line 
1#include "mltaln.h"
2
3#define DEBUG 0
4#define IODEBUG 0
5#define SCOREOUT 0
6
7#define END_OF_VEC -1
8
9int nadd;
10float thresholdtorev;
11int dodp;
12int addfragment;
13
14typedef struct _thread_arg
15{
16        int iend; 
17        char **seq;
18        char *tmpseq;
19        int *res;
20        int **spointt;
21        short *table1;
22        int iq;
23#ifdef enablemultithread
24        int *jshare;
25        int thread_no;
26        pthread_mutex_t *mutex_counter;
27#endif
28} thread_arg_t;
29
30
31
32void arguments( int argc, char *argv[] )
33{
34    int c;
35
36        nthread = 1;
37        inputfile = NULL;
38        nadd = 0;
39        dodp = 0;
40        alg = 'a';
41        alg = 'm';
42        dorp = NOTSPECIFIED;
43        fmodel = 0;
44//      ppenalty = (int)( -2.0 * 1000 - 0.5 );
45//      ppenalty_ex = (int)( -0.1 * 1000 - 0.5 );
46//      poffset = (int)( 0.1 * 1000 - 0.5 );
47        ppenalty = NOTSPECIFIED;
48        ppenalty_ex = NOTSPECIFIED;
49        poffset = NOTSPECIFIED;
50        kimuraR = 2;
51        pamN = 200;
52        thresholdtorev = 0.1;
53        addfragment = 0;
54
55
56    while( --argc > 0 && (*++argv)[0] == '-' )
57        {
58        while ( (c = *++argv[0]) )
59                {
60            switch( c )
61            {
62                                case 'i':
63                                        inputfile = *++argv;
64                                        fprintf( stderr, "inputfile = %s\n", inputfile );
65                                        --argc;
66                                        goto nextoption;
67                                case 'I':
68                                        nadd = atoi( *++argv );
69                                        fprintf( stderr, "nadd = %d\n", nadd );
70                                        --argc;
71                                        goto nextoption;
72                                case 'C':
73                                        nthread = atoi( *++argv );
74                                        fprintf( stderr, "nthread = %d\n", nthread );
75                                        --argc; 
76                                        goto nextoption;
77                                case 'f':
78                                        ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
79//                                      fprintf( stderr, "ppenalty = %d\n", ppenalty );
80                                        --argc;
81                                        goto nextoption;
82                                case 'g':
83                                        ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
84                                        fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
85                                        --argc;
86                                        goto nextoption;
87                                case 'h':
88                                        poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
89//                                      fprintf( stderr, "poffset = %d\n", poffset );
90                                        --argc;
91                                        goto nextoption;
92                                case 'k':
93                                        kimuraR = atoi( *++argv );
94                                        fprintf( stderr, "kappa = %d\n", kimuraR );
95                                        --argc;
96                                        goto nextoption;
97                                case 'j':
98                                        pamN = atoi( *++argv );
99                                        scoremtx = 0;
100                                        TMorJTT = JTT;
101                                        fprintf( stderr, "jtt/kimura %d\n", pamN );
102                                        --argc;
103                                        goto nextoption;
104                                case 't':
105                                        thresholdtorev = atof( *++argv );
106                                        fprintf( stderr, "thresholdtorev = %f\n", thresholdtorev );
107                                        --argc; 
108                                        goto nextoption;
109                                case 'd':
110                                        dodp = 1;
111                                        break;
112                                case 'F':
113                                        addfragment = 1;
114                                        break;
115#if 1
116                                case 'a':
117                                        fmodel = 1;
118                                        break;
119#endif
120                                case 'S':
121                                        alg = 'S';
122                                        break;
123                                case 'M':
124                                        alg = 'M';
125                                        break;
126                                case 'm':
127                                        alg = 'm';
128                                        break;
129                                case 'G':
130                                        alg = 'G';
131                                        break;
132                                case 'D':
133                                        dorp = 'd';
134                                        break;
135                                case 'P':
136                                        dorp = 'p';
137                                        break;
138                default:
139                    fprintf( stderr, "illegal option %c\n", c );
140                    argc = 0;
141                    break;
142            }
143                }
144                nextoption:
145                        ;
146        }
147    if( argc == 1 )
148    {
149        cut = atof( (*argv) );
150        argc--;
151    }
152    if( argc != 0 ) 
153    {
154        fprintf( stderr, "options: Check source file !\n" );
155        exit( 1 );
156    }
157        if( tbitr == 1 && outgap == 0 )
158        {
159                fprintf( stderr, "conflicting options : o, m or u\n" );
160                exit( 1 );
161        }
162}
163
164
165
166
167static int maxl;
168static int tsize;
169
170void seq_grp_nuc( int *grp, char *seq )
171{
172        int tmp;
173        int *grpbk = grp;
174        while( *seq )
175        {
176                tmp = amino_grp[(int)*seq++];
177                if( tmp < 4 )
178                        *grp++ = tmp;
179                else
180//                      fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
181                        ;
182        }
183        *grp = END_OF_VEC;
184        if( grp - grpbk < 6 )
185        {
186//              fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
187//              exit( 1 );
188                *grpbk = -1;
189        }
190}
191
192void seq_grp( int *grp, char *seq )
193{
194        int tmp;
195        int *grpbk = grp;
196        while( *seq )
197        {
198                tmp = amino_grp[(int)*seq++];
199                if( tmp < 6 )
200                        *grp++ = tmp;
201                else
202//                      fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
203                        ;
204        }
205        *grp = END_OF_VEC;
206        if( grp - grpbk < 6 )
207        {
208//              fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
209//              exit( 1 );
210                *grpbk = -1;
211        }
212}
213
214void makecompositiontable_p( short *table, int *pointt )
215{
216        int point;
217
218        while( ( point = *pointt++ ) != END_OF_VEC )
219                table[point]++;
220}
221
222
223void makepointtable_nuc( int *pointt, int *n )
224{
225        int point;
226        register int *p;
227
228        if( *n == -1 )
229        {
230                *pointt = -1;
231                return;
232        }
233
234        p = n;
235        point  = *n++ *  1024;
236        point += *n++ *   256;
237        point += *n++ *    64;
238        point += *n++ *    16;
239        point += *n++ *     4;
240        point += *n++;
241        *pointt++ = point;
242
243        while( *n != END_OF_VEC )
244        {
245                point -= *p++ * 1024;
246                point *= 4;
247                point += *n++;
248                *pointt++ = point;
249        }
250        *pointt = END_OF_VEC;
251}
252
253void makepointtable( int *pointt, int *n )
254{
255        int point;
256        register int *p;
257
258        if( *n == -1 )
259        {
260                *pointt = -1;
261                return;
262        }
263
264        p = n;
265        point  = *n++ *  7776;
266        point += *n++ *  1296;
267        point += *n++ *   216;
268        point += *n++ *    36;
269        point += *n++ *     6;
270        point += *n++;
271        *pointt++ = point;
272
273        while( *n != END_OF_VEC )
274        {
275                point -= *p++ * 7776;
276                point *= 6;
277                point += *n++;
278                *pointt++ = point;
279        }
280        *pointt = END_OF_VEC;
281}
282
283static int commonsextet_p2( short *table, int *pointt )
284{
285        int value = 0;
286        short tmp;
287        int point;
288        short *memo;
289        int *ct;
290        int *cp;
291
292        if( *pointt == -1 )
293                return( 0 );
294
295        memo = (short *)calloc( tsize, sizeof( short ) );
296        if( !memo ) ErrorExit( "Cannot allocate memo\n" );
297        ct = (int *)calloc( MIN( maxl, tsize )+1, sizeof( int ) ); // chuui!!
298        if( !ct ) ErrorExit( "Cannot allocate memo\n" );
299
300        cp = ct;
301        while( ( point = *pointt++ ) != END_OF_VEC )
302        {
303                tmp = memo[point]++;
304                if( tmp < table[point] )
305                        value++;
306                if( tmp == 0 ) *cp++ = point;
307        }
308        *cp = END_OF_VEC;
309       
310        cp =  ct;
311        while( *cp != END_OF_VEC )
312                memo[*cp++] = 0;
313
314        free( memo );
315        free( ct );
316        return( value );
317}
318
319static void     *directionthread( void *arg )
320{
321        thread_arg_t *targ = (thread_arg_t *)arg;
322        int iend = targ->iend;
323        char **seq = targ->seq;
324        char *tmpseq = targ->tmpseq;
325        int *res = targ->res;
326        int **spointt = targ->spointt;
327        short *table1 = targ->table1;
328        int iq = targ->iq;
329#ifdef enablemultithread
330        int thread_no = targ->thread_no;
331        int *jshare = targ->jshare; 
332#endif
333        int j;
334        char **mseq1, **mseq2;
335
336
337        if( dodp ) // nakuserukamo
338        {
339                mseq1 = AllocateCharMtx( 1, 0 );
340                mseq2 = AllocateCharMtx( 1, 0 );
341        }
342
343        j = -1;
344        while( 1 )
345        {
346#ifdef enablemultithread
347                if( nthread )
348                {
349                        pthread_mutex_lock( targ->mutex_counter );
350                        j = *jshare;
351                        if( j == iend )
352                        {
353                                fprintf( stderr, "\r %d / %d (thread %d)   \r", iq, njob, thread_no );
354                                pthread_mutex_unlock( targ->mutex_counter );
355                                break;
356                        }
357                        ++(*jshare);
358                        pthread_mutex_unlock( targ->mutex_counter );
359                }
360                else
361#endif
362                {
363                        j++;
364                        if( j == iend ) 
365                        {
366                                fprintf( stderr, "\r %d / %d  \r", iq, njob );
367                                break;
368                        }
369                }
370
371
372                if( dodp )
373                {
374//                      strcpy( mseq1[0], tmpseq );
375//                      strcpy( mseq2[0], seq[j] );
376                        mseq1[0] = tmpseq;
377                        mseq2[0] = seq[j];
378//                      res[j] = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, 0 );
379                        res[j] = L__align11_noalign( mseq1, mseq2 );
380                }
381                else
382                        res[j] = commonsextet_p2( table1, spointt[j] );
383        }
384        if( dodp ) // nakuserukamo
385        {
386                free( mseq1 );
387                free( mseq2 );
388//              G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
389                L__align11_noalign( NULL, NULL );
390        }
391//      else
392//              if( nthread )  // inthread == 0 no toki free suru to, error. nazeda
393//                      commonsextet_p( NULL, NULL );
394        return( NULL );
395}
396
397int main( int argc, char *argv[] )
398{
399        static int  *nlen;     
400        static int  *nogaplen; 
401        static char **name, **seq;
402        int i, j, istart, iend;
403        FILE *infp;
404//      FILE *adfp;
405        char c;
406
407        int *grpseq;
408        char *tmpseq, *revseq;
409        int  **pointt, **pointt_rev, **spointt;
410        float res_forward, res_reverse, res_max;
411        int ires, mres, mres2;
412        int *res;
413        static short *table1, *table1_rev;
414        static char **mseq1f, **mseq1r, **mseq2;
415
416        arguments( argc, argv );
417#ifndef enablemultithread
418        nthread = 0;
419#endif
420
421        if( inputfile )
422        {
423                infp = fopen( inputfile, "r" );
424                if( !infp )
425                {
426                        fprintf( stderr, "Cannot open %s\n", inputfile );
427                        exit( 1 );
428                }
429        }
430        else
431                infp = stdin;
432
433        getnumlen( infp );
434        rewind( infp );
435
436        if( alg == 'a' )
437        {
438                if( nlenmax < 10000 )
439                        alg = 'G';
440                else
441                        alg = 'S';
442        }
443
444        seq = AllocateCharMtx( njob, nlenmax*1+1 );
445
446#if 0
447        Read( name, nlen, seq );
448        readData( infp, name, nlen, seq );
449#else
450    name = AllocateCharMtx( njob, B+1 );
451    nlen = AllocateIntVec( njob ); 
452    nogaplen = AllocateIntVec( njob ); 
453        readData_pointer( infp, name, nlen, seq );
454        fclose( infp );
455
456        if( dorp != 'd' )
457        {
458                fprintf( stderr, "Not necessary!\n" );
459                for( i=0; i<njob; i++ ) 
460                        fprintf( stdout, "_F_%-10.10s\n", name[i]+1 );
461                exit( 1 );
462        }
463#endif
464
465        constants( njob, seq );
466
467
468#if 0
469        fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
470#endif
471
472        initSignalSM();
473
474        initFiles();
475
476        c = seqcheck( seq );
477        if( c )
478        {
479                fprintf( stderr, "Illegal character %c\n", c );
480                exit( 1 );
481        }
482
483        fprintf( stderr, "\n" );
484        if( alg == 'G' ) // do to the first sequence
485        {
486                mseq1f = AllocateCharMtx( 1, nlenmax+nlenmax );
487                mseq1r = AllocateCharMtx( 1, nlenmax+nlenmax );
488                mseq2 = AllocateCharMtx( 1, nlenmax+nlenmax );
489            tmpseq = AllocateCharVec( MAX( nlenmax, B ) +1 );
490
491                gappick0( mseq1f[0], seq[0] );
492                sreverse( mseq1r[0], mseq1f[0] );
493                strcpy( seq[0], mseq1f[0] );
494
495                if( nadd )
496                        istart = njob - nadd;
497                else
498                        istart = 1;
499
500                fprintf( stderr, "\n" );
501
502                for( i=0; i<istart; i++ )
503                {
504                        gappick0( tmpseq, seq[i] );
505                        strcpy( seq[i], tmpseq );
506                        strcpy( tmpseq, name[i] );
507                        strcpy( name[i], "_F_" );
508                        strncpy( name[i]+3, tmpseq+1, 10 );
509                        name[i][13] = 0;
510                }
511                for( i=istart; i<njob; i++ ) 
512                {
513                        gappick0( mseq2[0], seq[i] );
514
515//                      res_forward = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1f, mseq2, 0 );
516//                      res_reverse = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1r, mseq2, 0 );
517
518                        res_forward = L__align11_noalign( mseq1f, mseq2 );
519                        res_reverse = L__align11_noalign( mseq1r, mseq2 );
520#if 0
521
522                        strcpy( mseq2[0], seq[i] );
523                        strcpy( mseq1f[0], seq[0] );
524                        res_forward = G__align11( mseq1f, mseq2, nlenmax*2, 0, 0 );
525                        fprintf( stdout, "%s\n", mseq1f[0] );
526                        fprintf( stdout, "%s\n", mseq2[0] );
527
528                        strcpy( mseq2[0], seq[i] );
529                        sreverse( mseq1r[0], seq[0] );
530                        res_reverse = G__align11( mseq1r, mseq2, nlenmax*2, 0, 0 );
531                        fprintf( stdout, "%s\n", mseq1r[0] );
532                        fprintf( stdout, "%s\n", mseq2[0] );
533#endif
534
535//                      fprintf( stdout, "\nscore_for(%d,%d) = %f\n", 0, i, res_forward );
536//                      fprintf( stdout, "score_rev(%d,%d) = %f\n", 0, i, res_reverse );
537                        res_max = MAX(res_reverse,res_forward);
538                        if( (res_reverse-res_forward)/res_max > thresholdtorev ) // tekitou
539                        {
540//                              fprintf( stderr, "REVERSE!!!\n" );
541                                sreverse( seq[i], mseq2[0] );
542
543                                strcpy( tmpseq, name[i] );
544                                strcpy( name[i], "_R_" );
545                                strncpy( name[i]+3, tmpseq+1, 10 );
546                                name[i][13] = 0;
547                        }
548                        else
549                        {
550                                strcpy( seq[i], mseq2[0] );
551
552                                strcpy( tmpseq, name[i] );
553                                strcpy( name[i], "_F_" );
554                                strncpy( name[i]+3, tmpseq+1, 10 );
555                                name[i][13] = 0;
556                        }
557                }
558                FreeCharMtx( mseq1f );
559                FreeCharMtx( mseq1r );
560                FreeCharMtx( mseq2 );
561                free( tmpseq );
562        }
563        else if( alg == 'm' )
564        {
565                if( dodp ) // nakuserukamo
566                {
567                        mseq1f = AllocateCharMtx( 1, nlenmax+1);
568                        mseq1r = AllocateCharMtx( 1, nlenmax+1 );
569                        mseq2 = AllocateCharMtx( 1, nlenmax+1 );
570                }
571                else
572                {
573                        spointt = AllocateIntMtx( njob, 0 ); 
574                        pointt = AllocateIntMtx( njob, nlenmax+1 );
575                        pointt_rev = AllocateIntMtx( njob, nlenmax+1 );
576                }
577            tmpseq = AllocateCharVec( MAX( nlenmax, B ) +1 );
578            revseq = AllocateCharVec( nlenmax+1 );
579                grpseq = AllocateIntVec( nlenmax+1 );
580                res = AllocateIntVec( njob );
581                if( dorp == 'd' ) tsize = (int)pow( 4, 6 );
582                else              tsize = (int)pow( 6, 6 ); // iranai
583
584                maxl = 0;
585                for( i=0; i<njob; i++ )
586                {
587                        gappick0( tmpseq, seq[i] );
588                        nogaplen[i] = strlen( tmpseq );
589                        if( nogaplen[i] > maxl ) maxl = nogaplen[i];
590                }
591
592                if( nadd )
593                        iend = njob - nadd;
594                else
595                        iend = 1;
596
597                for( i=0; i<iend; i++ )
598                {
599//                      fprintf( stdout, "%d, SKIP\n", i );
600                        gappick0( tmpseq, seq[i] );
601                        strcpy( seq[i], tmpseq );
602//                      if( !nadd ) strcpy( seq[i], tmpseq ); // seq ha tsukawanaikara ii.
603
604                        if( !dodp )
605                        {
606                                seq_grp_nuc( grpseq, tmpseq );
607                                makepointtable_nuc( pointt[i], grpseq );
608                                spointt[i] = pointt[i];
609                        }
610
611                        strcpy( tmpseq, name[i] );
612                        strcpy( name[i], "_F_" );
613                        strncpy( name[i]+3, tmpseq+1, 10 );
614                        name[i][13] = 0;
615                }
616
617                if( nadd )
618                        istart = njob - nadd;
619                else
620                        istart = 1;
621
622                fprintf( stderr, "\n" );
623                for( i=istart; i<njob; i++ ) 
624                {
625//                      fprintf( stderr, "\r %d / %d ", i, njob );
626                        gappick0( tmpseq, seq[i] );
627                        strcpy( seq[i], tmpseq );
628                        sreverse( revseq, tmpseq );
629
630                        if( !dodp )
631                        {
632                                table1 = (short *)calloc( tsize, sizeof( short ) );
633                                if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
634                                table1_rev = (short *)calloc( tsize, sizeof( short ) );
635                                if( !table1_rev ) ErrorExit( "Cannot allocate table1_rev\n" );
636                                seq_grp_nuc( grpseq, tmpseq );
637                                makepointtable_nuc( pointt[i], grpseq );
638                                makecompositiontable_p( table1, pointt[i] );
639                                seq_grp_nuc( grpseq, revseq );
640                                makepointtable_nuc( pointt_rev[i], grpseq );
641                                makecompositiontable_p( table1_rev, pointt_rev[i] );
642                        }
643
644                        if( nadd && addfragment )
645                                iend = njob-nadd;
646                        else
647                                iend = i;
648
649
650#ifdef enablemultithread
651                        if( nthread )
652                        {
653                                pthread_t *handle;
654                                pthread_mutex_t mutex_counter;
655                                thread_arg_t *targ;
656                                int *jsharept;
657               
658                                targ = calloc( nthread, sizeof( thread_arg_t ) );
659                                handle = calloc( nthread, sizeof( pthread_t ) );
660                                pthread_mutex_init( &mutex_counter, NULL );
661                                jsharept = calloc( 1, sizeof(int) );
662                                *jsharept = 0;
663               
664                                for( j=0; j<nthread; j++ )
665                                {
666                                        targ[j].iend = iend;
667                                        targ[j].seq = seq;
668                                        targ[j].tmpseq = tmpseq; 
669                                        targ[j].res = res; 
670                                        targ[j].spointt = spointt; 
671                                        targ[j].table1 = table1; 
672                                        targ[j].jshare = jsharept;
673                                        targ[j].iq = i;
674                                        targ[j].mutex_counter = &mutex_counter;
675                                        targ[j].thread_no = j;
676                                        pthread_create( handle+j, NULL, directionthread, (void *)(targ+j) );
677                                }
678                                for( j=0; j<nthread; j++ ) pthread_join( handle[j], NULL );
679                                pthread_mutex_destroy( &mutex_counter );
680                                free( handle );
681                                free( targ );
682                                free( jsharept );
683                        }
684                        else
685#endif
686                        {
687                                thread_arg_t *targ;
688                                targ = calloc( 1, sizeof( thread_arg_t ) );
689                                targ[0].iend = iend;
690                                targ[0].seq = seq;
691                                targ[0].tmpseq = tmpseq; 
692                                targ[0].res = res; 
693                                targ[0].spointt = spointt; 
694                                targ[0].table1 = table1; 
695                                targ[0].iq = i; 
696                                directionthread( targ );
697                                free( targ );
698                        }
699
700
701                        mres = mres2 = 0;
702                        for( j=0; j<iend; j++ )
703                        {
704                                ires = res[j];
705//                              fprintf( stdout, "ires (%d,%d) = %d\n", i, j, ires );
706//                              fflush( stdout );
707                                if( ires>mres2 ) 
708                                {
709                                        if( ires>mres ) 
710                                        {
711                                                mres2 = mres;
712                                                mres = ires;
713                                        }
714                                        else
715                                                mres2 = ires;
716                                }
717                        }
718                        res_forward = (float)( mres + mres2 ) / 2;
719
720#ifdef enablemultithread
721                        if( nthread )
722                        {
723                                pthread_t *handle;
724                                pthread_mutex_t mutex_counter;
725                                thread_arg_t *targ;
726                                int *jsharept;
727               
728                                targ = calloc( nthread, sizeof( thread_arg_t ) );
729                                handle = calloc( nthread, sizeof( pthread_t ) );
730                                pthread_mutex_init( &mutex_counter, NULL );
731                                jsharept = calloc( 1, sizeof(int) );
732                                *jsharept = 0;
733               
734                                for( j=0; j<nthread; j++ )
735                                {
736                                        targ[j].iend = iend;
737                                        targ[j].seq = seq;
738                                        targ[j].tmpseq = revseq; 
739                                        targ[j].res = res; 
740                                        targ[j].spointt = spointt; 
741                                        targ[j].table1 = table1_rev; 
742                                        targ[j].jshare = jsharept;
743                                        targ[j].iq = i;
744                                        targ[j].mutex_counter = &mutex_counter;
745                                        targ[j].thread_no = j;
746                                        pthread_create( handle+j, NULL, directionthread, (void *)(targ+j) );
747                                }
748                                for( j=0; j<nthread; j++ ) pthread_join( handle[j], NULL );
749                                pthread_mutex_destroy( &mutex_counter );
750                                free( handle );
751                                free( targ );
752                                free( jsharept );
753                        }
754                        else
755#endif
756                        {
757                                thread_arg_t *targ;
758                                targ = calloc( 1, sizeof( thread_arg_t ) );
759                                targ[0].iend = iend;
760                                targ[0].seq = seq;
761                                targ[0].tmpseq = revseq; 
762                                targ[0].res = res; 
763                                targ[0].spointt = spointt;
764                                targ[0].table1 = table1_rev; 
765                                targ[0].iq = i; 
766                                directionthread( targ );
767                                free( targ );
768                        }
769
770                        mres = mres2 = 0;
771                        for( j=0; j<iend; j++ )
772                        {
773                                ires = res[j];
774                                if( ires>mres2 )
775                                {
776                                        if( ires>mres ) 
777                                        {
778                                                mres2 = mres;
779                                                mres = ires;
780                                        }
781                                        else
782                                                mres2 = ires;
783                                }
784                        }
785                        res_reverse = (float)( mres + mres2 ) / 2;
786
787//                      fprintf( stdout, "\n" );
788//                      fprintf( stdout, "score_for(%d,%d) = %f\n", 0, i, res_forward );
789//                      fprintf( stdout, "score_rev(%d,%d) = %f\n", 0, i, res_reverse );
790//                      fflush( stdout );
791                        res_max = MAX(res_reverse,res_forward);
792                        if( (res_reverse-res_forward)/res_max > thresholdtorev ) // tekitou
793
794                        {
795                                strcpy( seq[i], revseq );
796
797                                strcpy( tmpseq, name[i] );
798                                strcpy( name[i], "_R_" );
799                                strncpy( name[i]+3, tmpseq+1, 10 );
800                                name[i][13] = 0;
801                                if( !dodp ) spointt[i] = pointt_rev[i];
802                        }
803                        else
804                        {
805                                strcpy( tmpseq, name[i] );
806                                strcpy( name[i], "_F_" );
807                                strncpy( name[i]+3, tmpseq+1, 10 );
808                                name[i][13] = 0;
809                                if( !dodp ) spointt[i] = pointt[i];
810                        }
811
812                        if( !dodp )
813                        {
814                                free( table1 );
815                                free( table1_rev );
816                        }
817                }
818
819                free( grpseq );
820                free( tmpseq );
821                free( revseq );
822                free( res );
823                if( dodp )
824                {
825                        FreeCharMtx( mseq1f );
826                        FreeCharMtx( mseq1r );
827                        FreeCharMtx( mseq2 );
828                }
829                else
830                {
831                        FreeIntMtx( pointt );
832                        FreeIntMtx( pointt_rev );
833                        free( spointt );
834                }
835        }
836        else
837        {
838                fprintf( stderr, "Unknown alg %c\n", alg );
839                exit( 1 );
840        }
841//      writeData_pointer( stdout, njob, name, nlen, seq );
842        for( i=0; i<njob; i++ ) 
843        {
844//              fprintf( stdout, ">%s\n", name[i] );
845//              fprintf( stdout, "%s\n", seq[i] );
846                fprintf( stdout, "%s\n", name[i] );
847        }
848
849        fprintf( stderr, "\n" );
850        SHOWVERSION;
851        return( 0 );
852}
853
Note: See TracBrowser for help on using the repository browser.