source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/dndpre2.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: 3.6 KB
Line 
1#include "mltaln.h"
2
3#define DEBUG 0
4#define IODEBUG 0
5#define SCOREOUT 1
6
7
8void arguments( int argc, char *argv[] )
9{
10    int c;
11
12        nblosum = 62;
13        fmodel = 0;
14        calledByXced = 0;
15        scoremtx = 0;
16        dorp = NOTSPECIFIED;
17        ppenalty = NOTSPECIFIED;
18        ppenalty_ex = NOTSPECIFIED;
19        poffset = NOTSPECIFIED;
20        kimuraR = NOTSPECIFIED;
21        pamN = NOTSPECIFIED;
22        geta2 = GETA2;
23        fftWinSize = NOTSPECIFIED;
24        fftThreshold = NOTSPECIFIED;
25
26    while( --argc > 0 && (*++argv)[0] == '-' )
27        {
28        while ( (c = *++argv[0]) )
29                {
30            switch( c )
31            {
32                                case 'f':
33                                        ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
34                                        fprintf( stderr, "ppenalty = %d\n", ppenalty );
35                                        --argc;
36                                        goto nextoption;
37                                case 'g':
38                                        ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
39                                        fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
40                                        --argc;
41                                        goto nextoption;
42                                case 'h':
43                                        poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
44                                        fprintf( stderr, "poffset = %d\n", poffset );
45                                        --argc;
46                                        goto nextoption;
47                                case 'k':
48                                        kimuraR = atoi( *++argv );
49                                        fprintf( stderr, "kimuraR = %d\n", kimuraR );
50                                        --argc;
51                                        goto nextoption;
52                                case 'b':
53                                        nblosum = atoi( *++argv );
54                                        scoremtx = 1;
55                                        fprintf( stderr, "blosum %d\n", nblosum );
56                                        --argc;
57                                        goto nextoption;
58                                case 'j':
59                                        pamN = atoi( *++argv );
60                                        scoremtx = 0;
61                                        fprintf( stderr, "jtt %d\n", pamN );
62                                        --argc;
63                                        goto nextoption;
64                                case 'm':
65                                        fmodel = 1;
66                                        break;
67                                case 'r':
68                                        fmodel = -1;
69                                        break;
70                                case 'D':
71                                        dorp = 'd';
72                                        break;
73                                case 'P':
74                                        dorp = 'p';
75                                        break;
76                default:
77                    fprintf( stderr, "illegal option %c\n", c );
78                    argc = 0;
79                    break;
80            }
81                }
82                nextoption:
83                        ;
84        }
85    if( argc == 1 )
86    {
87        cut = atof( (*argv) );
88        argc--;
89    }
90    if( argc != 0 ) 
91    {
92        fprintf( stderr, "options: Check source file !\n" );
93        exit( 1 );
94    }
95        if( tbitr == 1 && outgap == 0 )
96        {
97                fprintf( stderr, "conflicting options : o, m or u\n" );
98                exit( 1 );
99        }
100        if( alg == 'C' && outgap == 0 )
101        {
102                fprintf( stderr, "conflicting options : C, o\n" );
103                exit( 1 );
104        }
105}
106
107int main( int argc, char *argv[] )
108{
109        static int  nlen[M];   
110        static char **name, **seq;
111        int i, j, alloclen, c;
112        double **mtx;
113        double *self;
114        double tmpdouble;
115        FILE *fp;
116
117        arguments( argc, argv );
118
119        getnumlen( stdin );
120        rewind( stdin );
121
122        if( njob < 2 )
123        {
124                fprintf( stderr, "At least 2 sequences should be input!\n"
125                                                 "Only %d sequence found.\n", njob ); 
126                exit( 1 );
127        }
128
129        name = AllocateCharMtx( njob, B+1 );
130        seq = AllocateCharMtx( njob, nlenmax*9+1 );
131        mtx = AllocateDoubleMtx( njob, njob );
132        self = AllocateDoubleVec( njob );
133        alloclen = nlenmax*9;
134
135        readData_pointer( stdin, name, nlen, seq );
136        constants( njob, seq );
137
138
139
140
141        c = seqcheck( seq );
142        if( c )
143        {
144                fprintf( stderr, "Illeagal character %c\n", c );
145                exit( 1 );
146        }
147
148        for( i=0; i<njob; i++ ) 
149        {
150                self[i] = (double)substitution_nid( seq[i], seq[i] );
151//              fprintf( stdout, "self[%d] = %f\n", i, self[i] );
152        }
153
154        for( i=0; i<njob-1; i++ ) 
155                for( j=i+1; j<njob; j++ ) 
156                {
157                        tmpdouble = (double)substitution_score( seq[i], seq[j] );
158//                      fprintf( stdout, "tmpdouble = %f\n", tmpdouble );
159                        mtx[i][j] = ( 1.0 - tmpdouble / MIN( self[i], self[j] ) );
160                        if( mtx[i][j] < 0.95 )
161                                mtx[i][j] = - log( 1.0 - mtx[i][j] );
162                        else
163                                mtx[i][j] = 3.0;
164                }
165       
166#if TEST
167        for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ ) 
168                fprintf( stdout, "i=%d, j=%d, mtx[][] = %f\n", i, j, mtx[i][j] );
169#endif
170
171        fp = fopen( "hat2", "w" );
172        WriteHat2_pointer( fp, njob, name, mtx );
173        fclose( fp );
174        exit( 0 );
175
176        return( 0 );
177}
Note: See TracBrowser for help on using the repository browser.