source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/sextet5.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: 5.9 KB
Line 
1#include "mltaln.h"
2#include "mtxutl.h"
3
4#define DEBUG 0
5#define TEST  0
6
7#define END_OF_VEC -1
8
9static int maxl;
10static int tsize;
11
12void arguments( int argc, char *argv[] )
13{
14        int c;
15
16        inputfile = NULL;
17        disopt = 0;
18        scoremtx = 1;
19        nblosum = 62;
20        dorp = NOTSPECIFIED;
21
22    while( --argc > 0 && (*++argv)[0] == '-' )
23        {
24        while ( ( c = *++argv[0] ) )
25                {
26            switch( c )
27            {
28                                case 'i':
29                                        inputfile = *++argv;
30                                        fprintf( stderr, "inputfile = %s\n", inputfile );
31                                        --argc;
32                                        goto nextoption;
33                                case 'D':
34                                        dorp = 'd';
35                                        break;
36                                case 'P':
37                                        dorp = 'p';
38                                        break;
39                                case 'I':
40                                        disopt = 1;
41                                        break;
42                default:
43                    fprintf( stderr, "illegal option %c\n", c );
44                    argc = 0;
45                    break;
46            }
47                }
48                nextoption:
49                        ;
50        }
51    if( argc != 0 )
52    {
53        fprintf( stderr, "options: -i\n" );
54        exit( 1 );
55    }
56}
57
58void seq_grp_nuc( int *grp, char *seq )
59{
60        int tmp;
61        while( *seq )
62        {
63                tmp = amino_grp[(int)*seq++];
64                if( tmp < 4 )
65                        *grp++ = tmp;
66                else
67                        fprintf( stderr, "WARNING : Unknown character %c\n", *(seq-1) );
68        }
69        *grp = END_OF_VEC;
70}
71
72void seq_grp( int *grp, char *seq )
73{
74        int tmp;
75        while( *seq )
76        {
77                tmp = amino_grp[(int)*seq++];
78                if( tmp < 6 )
79                        *grp++ = tmp;
80                else
81                        fprintf( stderr, "WARNING : Unknown character %c\n", *(seq-1) );
82        }
83        *grp = END_OF_VEC;
84}
85
86void makecompositiontable_p( short *table, int *pointt )
87{
88        int point;
89
90        while( ( point = *pointt++ ) != END_OF_VEC )
91                table[point]++;
92}
93
94int commonsextet_p( short *table, int *pointt )
95{
96        int value = 0;
97        short tmp;
98        int point;
99        static short *memo = NULL;
100        static int *ct = NULL;
101        static int *cp;
102
103        if( !memo )
104        {
105                memo = (short *)calloc( tsize, sizeof( short ) );
106                if( !memo ) ErrorExit( "Cannot allocate memo\n" );
107                ct = (int *)calloc( MIN( maxl, tsize)+1, sizeof( int ) );
108                if( !ct ) ErrorExit( "Cannot allocate memo\n" );
109        }
110
111        cp = ct;
112        while( ( point = *pointt++ ) != END_OF_VEC )
113        {
114                tmp = memo[point]++;
115                if( tmp < table[point] )
116                        value++;
117                if( tmp == 0 ) *cp++ = point;
118//              fprintf( stderr, "cp - ct = %d (tsize = %d)\n", cp - ct, tsize );
119        }
120        *cp = END_OF_VEC;
121       
122        cp =  ct;
123        while( *cp != END_OF_VEC )
124                memo[*cp++] = 0;
125
126        return( value );
127}
128
129void makepointtable_nuc( int *pointt, int *n )
130{
131        int point;
132        register int *p;
133
134        p = n;
135        point  = *n++ *  1024;
136        point += *n++ *   256;
137        point += *n++ *    64;
138        point += *n++ *    16;
139        point += *n++ *     4;
140        point += *n++;
141        *pointt++ = point;
142
143        while( *n != END_OF_VEC )
144        {
145                point -= *p++ * 1024;
146                point *= 4;
147                point += *n++;
148                *pointt++ = point;
149        }
150        *pointt = END_OF_VEC;
151}
152
153void makepointtable( int *pointt, int *n )
154{
155        int point;
156        register int *p;
157
158        p = n;
159        point  = *n++ *  7776;
160        point += *n++ *  1296;
161        point += *n++ *   216;
162        point += *n++ *    36;
163        point += *n++ *     6;
164        point += *n++;
165        *pointt++ = point;
166
167        while( *n != END_OF_VEC )
168        {
169                point -= *p++ * 7776;
170                point *= 6;
171                point += *n++;
172                *pointt++ = point;
173        }
174        *pointt = END_OF_VEC;
175}
176
177int main( int argc, char **argv )
178{
179        int i, j;
180        FILE *fp, *infp;
181        char **seq;
182        int *grpseq;
183        char *tmpseq;
184        int  **pointt;
185        static char **name;
186        static int nlen[M];
187        double **mtx;
188        double **mtx2;
189        double score, score0;
190        static short *table1;
191        char b[B];
192
193        arguments( argc, argv );
194
195        if( inputfile )
196        {
197                infp = fopen( inputfile, "r" );
198                if( !infp )
199                {
200                        fprintf( stderr, "Cannot open %s\n", inputfile );
201                        exit( 1 );
202                }
203        }
204        else
205                infp = stdin;
206
207#if 0
208        PreRead( stdin, &njob, &nlenmax );
209#else
210        getnumlen( infp );
211#endif
212        rewind( infp );
213        if( njob < 2 )
214        {
215                fprintf( stderr, "At least 2 sequences should be input!\n"
216                                                 "Only %d sequence found.\n", njob );
217                exit( 1 );
218        }
219
220        name = AllocateCharMtx( njob, B+1 );
221        tmpseq = AllocateCharVec( nlenmax+1 );
222        seq = AllocateCharMtx( njob, nlenmax+1 );
223        grpseq = AllocateIntVec( nlenmax+1 );
224        pointt = AllocateIntMtx( njob, nlenmax+1 );
225        mtx = AllocateDoubleMtx( njob, njob );
226        mtx2 = AllocateDoubleMtx( njob, njob );
227        pamN = NOTSPECIFIED;
228
229#if 0
230        FRead( infp, name, nlen, seq );
231#else
232        readData_pointer( infp, name, nlen, seq );
233#endif
234
235        fclose( infp );
236
237        constants( njob, seq );
238
239        if( dorp == 'd' ) tsize = (int)pow( 4, 6 );
240        else              tsize = (int)pow( 6, 6 );
241
242        maxl = 0;
243        for( i=0; i<njob; i++ ) 
244        {
245                gappick0( tmpseq, seq[i] );
246                nlen[i] = strlen( tmpseq );
247                if( nlen[i] < 6 )
248                {
249                        fprintf( stderr, "Seq %d, too short, %d characters\n", i+1, nlen[i] );
250                        exit( 1 );
251                }
252                if( nlen[i] > maxl ) maxl = nlen[i];
253                if( dorp == 'd' ) /* nuc */
254                {
255                        seq_grp_nuc( grpseq, tmpseq );
256                        makepointtable_nuc( pointt[i], grpseq );
257                }
258                else                 /* amino */
259                {
260                        seq_grp( grpseq, tmpseq );
261                        makepointtable( pointt[i], grpseq );
262                }
263        }
264        for( i=0; i<njob; i++ )
265        {
266                table1 = (short *)calloc( tsize, sizeof( short ) );
267                if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
268                if( i % 10 == 0 )
269                {
270                        fprintf( stderr, "%4d / %4d\r", i+1, njob );
271                }
272                makecompositiontable_p( table1, pointt[i] );
273
274                for( j=i; j<njob; j++ ) 
275                {
276                        score = (double)commonsextet_p( table1, pointt[j] );
277                        mtx[i][j] = score;
278                } 
279                free( table1 );
280        }
281        for( i=0; i<njob; i++ )
282        {
283                score0 = mtx[i][i];
284                for( j=0; j<njob; j++ ) 
285                        mtx2[i][j] = ( score0 - mtx[MIN(i,j)][MAX(i,j)] ) / score0 * 3.0;
286        }
287        for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ ) 
288        {
289#if TEST
290                double jscore;
291                jscore = mtx[i][j] / ( MIN( strlen( seq[i] ), strlen( seq[j] ) ) - 2 );
292                fprintf( stdout, "jscore = %f\n", jscore );
293
294                fprintf( stdout, "mtx2[%d][%d] = %f, mtx2[%d][%d] = %f\n", i, j, mtx2[i][j], j, i, mtx2[j][i] );
295#endif
296                mtx2[i][j] = MIN( mtx2[i][j], mtx2[j][i] );
297#if TEST
298                fprintf( stdout, "sonokekka mtx2[%d][%d] %f\n", i, j, mtx2[i][j] );
299#endif
300        }
301
302        if( disopt )
303        {
304                for( i=0; i<njob; i++ ) 
305                {
306                        sprintf( b, "=lgth = %04d", nlen[i] );
307                        strins( b, name[i] );
308                }
309        }
310               
311        fp = fopen( "hat2", "w" );
312        if( !fp ) ErrorExit( "Cannot open hat2." );
313        WriteHat2_pointer( fp, njob, name, mtx2 );
314        fclose( fp );
315
316        fprintf( stderr, "\n" );
317        SHOWVERSION;
318        exit( 0 );
319}
Note: See TracBrowser for help on using the repository browser.