5 static char *whereiscontrafold;
7 void unknown_n( char *out, char *in )
11 if( *in == 'a' || *in == 'A' )
13 else if( *in == 't' || *in == 'T' || *in == 'u' || *in == 'U' )
15 else if( *in == 'g' || *in == 'G' )
17 else if( *in == 'c' || *in == 'C' )
30 void outcontrafold( FILE *fp, RNApair **pairprob, int length )
34 for( i=0; i<length; i++ ) for( pt=pairprob[i]; pt->bestpos!=-1; pt++ )
37 fprintf( fp, "%d %d %f\n", i, pt->bestpos, pt->bestscore );
42 static void readcontrafold( FILE *fp, RNApair **pairprob, int length )
51 pairnum = (int *)calloc( length, sizeof( int ) );
52 for( i=0; i<length; i++ ) pairnum[i] = 0;
56 if( feof( fp ) ) break;
57 fgets( gett, 9999, fp );
59 // fprintf( stderr, "gett=%s\n", gett );
63 sscanf( gett, "%d ", &left );
66 // fprintf( stderr, "left=%d\n", left );
67 pt = strchr( pt, ' ' ) + 1;
68 // fprintf( stderr, "pt=%s\n", pt );
70 while( (pt = strchr( pt, ' ' ) ) )
73 // fprintf( stderr, "pt=%s\n", pt );
74 sscanf( pt, "%d:%f", &right, &prob );
77 // fprintf( stderr, "%d-%d, %f\n", left, right, prob );
79 pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
80 pairprob[left][pairnum[left]].bestscore = prob;
81 pairprob[left][pairnum[left]].bestpos = right;
83 pairprob[left][pairnum[left]].bestscore = -1.0;
84 pairprob[left][pairnum[left]].bestpos = -1;
85 // fprintf( stderr, "%d-%d, %f\n", left, right, prob );
87 pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
88 pairprob[right][pairnum[right]].bestscore = prob;
89 pairprob[right][pairnum[right]].bestpos = left;
91 pairprob[right][pairnum[right]].bestscore = -1.0;
92 pairprob[right][pairnum[right]].bestpos = -1;
93 // fprintf( stderr, "%d-%d, %f\n", right, left, prob );
100 void arguments( int argc, char *argv[] )
105 kimuraR = NOTSPECIFIED;
107 whereiscontrafold = NULL;
109 while( --argc > 0 && (*++argv)[0] == '-' )
111 while ( (c = *++argv[0]) )
117 fprintf( stderr, "inputfile = %s\n", inputfile );
121 whereiscontrafold = *++argv;
122 fprintf( stderr, "whereiscontrafold = %s\n", whereiscontrafold );
126 fprintf( stderr, "illegal option %c\n", c );
136 fprintf( stderr, "options: Check source file !\n" );
142 int main( int argc, char *argv[] )
144 static char com[10000];
148 static char **name, **seq, **nogap;
154 RNApair **alnpairprob;
161 arguments( argc, argv );
165 infp = fopen( inputfile, "r" );
168 fprintf( stderr, "Cannot open %s\n", inputfile );
175 if( !whereiscontrafold )
176 whereiscontrafold = "";
183 fprintf( stderr, "nuc only\n" );
187 seq = AllocateCharMtx( njob, nlenmax*2+1 );
188 nogap = AllocateCharMtx( njob, nlenmax*2+1 );
189 gapmap = AllocateIntMtx( njob, nlenmax*2+1 );
190 order = AllocateIntVec( njob );
191 name = AllocateCharMtx( njob, B+1 );
192 nlen = AllocateIntVec( njob );
193 pairprob = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
194 alnpairprob = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
195 alnpairnum = AllocateIntVec( nlenmax );
197 for( i=0; i<nlenmax; i++ ) alnpairnum[i] = 0;
199 readData_pointer( infp, name, nlen, seq );
202 for( i=0; i<njob; i++ )
204 pairprob[i] = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
205 for( j=0; j<nlenmax; j++ )
207 pairprob[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
208 pairprob[i][j][0].bestpos = -1;
209 pairprob[i][j][0].bestscore = -1.0;
211 unknown_n( nogap[i], seq[i] );
214 for( j=0; j<nlenmax; j++ )
216 alnpairprob[j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
217 alnpairprob[j][0].bestpos = -1;
218 alnpairprob[j][0].bestscore = -1.0;
222 constants( njob, seq );
224 fprintf( stderr, "running contrafold\n" );
225 for( i=0; i<njob; i++ )
227 fprintf( stderr, "%d / %d\n", i+1, njob );
228 commongappick_record( 1, nogap+i, gapmap[i] );
229 infp = fopen( "_contrafoldin", "w" );
230 fprintf( infp, ">in\n%s\n", nogap[i] );
232 #if 0 // contrafold v1
233 sprintf( com, "env PATH=%s contrafold predict _contrafoldin --posteriors 0.01 > _contrafoldout", whereiscontrafold );
234 #else // contrafold v2
235 sprintf( com, "env PATH=%s contrafold predict _contrafoldin --posteriors 0.01 _contrafoldout", whereiscontrafold );
240 fprintf( stderr, "error in contrafold\n" );
241 fprintf( stderr, "=================================================================\n" );
242 fprintf( stderr, "=================================================================\n" );
243 fprintf( stderr, "==\n" );
244 fprintf( stderr, "== This version of MAFFT supports CONTRAfold v2.02.\n" );
245 fprintf( stderr, "== If you have a lower version of CONTRAfold installed in the\n" );
246 fprintf( stderr, "== %s directory,\n", whereiscontrafold );
247 fprintf( stderr, "== please update it!\n" );
248 fprintf( stderr, "==\n" );
249 fprintf( stderr, "=================================================================\n" );
250 fprintf( stderr, "=================================================================\n" );
255 infp = fopen( "_contrafoldout", "r" );
256 readcontrafold( infp, pairprob[i], nlenmax );
258 fprintf( stdout, ">%d\n", i );
259 outcontrafold( stdout, pairprob[i], nlenmax );
262 for( i=0; i<njob; i++ )
264 for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
267 right = gapmap[i][pairprobpt->bestpos];
268 prob = pairprobpt->bestscore;
270 for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
271 if( pt->bestpos == right ) break;
273 if( pt->bestpos == -1 )
275 alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
276 adpos = alnpairnum[left];
278 alnpairprob[left][adpos].bestscore = 0.0;
279 alnpairprob[left][adpos].bestpos = right;
280 alnpairprob[left][adpos+1].bestscore = -1.0;
281 alnpairprob[left][adpos+1].bestpos = -1;
282 pt = alnpairprob[left]+adpos;
285 adpos = pt-alnpairprob[left];
287 pt->bestscore += prob;
288 if( pt->bestpos != right )
290 fprintf( stderr, "okashii!\n" );
293 // fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
299 fprintf( stdout, "result=\n" );
301 for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
303 pairprobpt->bestscore /= (float)njob;
305 right = pairprobpt->bestpos;
306 prob = pairprobpt->bestscore;
307 fprintf( stdout, "%d-%d, %f\n", left, right, prob );