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 sprintf( com, "env PATH=%s contrafold predict _contrafoldin --posteriors 0.01 > _contrafoldout", whereiscontrafold );
236 fprintf( stderr, "error in contrafold\n" );
241 infp = fopen( "_contrafoldout", "r" );
242 readcontrafold( infp, pairprob[i], nlenmax );
244 fprintf( stdout, ">%d\n", i );
245 outcontrafold( stdout, pairprob[i], nlenmax );
248 for( i=0; i<njob; i++ )
250 for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
253 right = gapmap[i][pairprobpt->bestpos];
254 prob = pairprobpt->bestscore;
256 for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
257 if( pt->bestpos == right ) break;
259 if( pt->bestpos == -1 )
261 alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
262 adpos = alnpairnum[left];
264 alnpairprob[left][adpos].bestscore = 0.0;
265 alnpairprob[left][adpos].bestpos = right;
266 alnpairprob[left][adpos+1].bestscore = -1.0;
267 alnpairprob[left][adpos+1].bestpos = -1;
268 pt = alnpairprob[left]+adpos;
271 adpos = pt-alnpairprob[left];
273 pt->bestscore += prob;
274 if( pt->bestpos != right )
276 fprintf( stderr, "okashii!\n" );
279 // fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
285 fprintf( stdout, "result=\n" );
287 for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
289 pairprobpt->bestscore /= (float)njob;
291 right = pairprobpt->bestpos;
292 prob = pairprobpt->bestscore;
293 fprintf( stdout, "%d-%d, %f\n", left, right, prob );