source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/rna.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: 15.1 KB
Line 
1#include "mltaln.h"
2#include "dp.h"
3
4#define MEMSAVE 1
5
6#define DEBUG 1
7#define USE_PENALTY_EX  1
8#define STOREWM 1
9
10
11
12static float singleribosumscore( int n1, int n2, char **s1, char **s2, double *eff1, double *eff2, int p1, int p2 )
13{
14        float val;
15        int i, j;
16        int code1, code2;
17
18        val = 0.0;
19        for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
20        {
21                code1 = amino_n[(int)s1[i][p1]];
22                if( code1 > 3 ) code1 = 36;
23                code2 = amino_n[(int)s2[j][p2]];
24                if( code2 > 3 ) code2 = 36;
25
26//              fprintf( stderr, "'l'%c-%c: %f\n", s1[i][p1], s2[j][p2], (float)ribosumdis[code1][code2] );
27
28                val += (float)ribosumdis[code1][code2] * eff1[i] * eff2[j];
29        }
30        return( val );
31}
32static float pairedribosumscore53( int n1, int n2, char **s1, char **s2, double *eff1, double *eff2, int p1, int p2, int c1, int c2 )
33{
34        float val;
35        int i, j;
36        int code1o, code1u, code2o, code2u, code1, code2;
37
38        val = 0.0;
39        for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
40        {
41                code1o = amino_n[(int)s1[i][p1]];
42                code1u = amino_n[(int)s1[i][c1]];
43                if( code1o > 3 ) code1 = code1o = 36;
44                else if( code1u > 3 ) code1 = 36;
45                else code1 = 4 + code1o * 4 + code1u;
46
47                code2o = amino_n[(int)s2[j][p2]];
48                code2u = amino_n[(int)s2[j][c2]];
49                if( code2o > 3 ) code2 = code1o = 36;
50                else if( code2u > 3 ) code2 = 36;
51                else code2 = 4 + code2o * 4 + code2u;
52
53
54//              fprintf( stderr, "%c%c-%c%c: %f\n", s1[i][p1], s1[i][c1], s2[j][p2], s2[j][c2], (float)ribosumdis[code1][code2] );
55
56                if( code1 == 36 || code2 == 36 )
57                        val += (float)n_dis[code1o][code2o] * eff1[i] * eff2[j];
58                else
59                        val += (float)ribosumdis[code1][code2] * eff1[i] * eff2[j];
60        }
61        return( val );
62}
63
64static float pairedribosumscore35( int n1, int n2, char **s1, char **s2, double *eff1, double *eff2, int p1, int p2, int c1, int c2 )
65{
66        float val;
67        int i, j;
68        int code1o, code1u, code2o, code2u, code1, code2;
69
70        val = 0.0;
71        for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
72        {
73                code1o = amino_n[(int)s1[i][p1]];
74                code1u = amino_n[(int)s1[i][c1]];
75                if( code1o > 3 ) code1 = code1o = 36;
76                else if( code1u > 3 ) code1 = 36;
77                else code1 = 4 + code1u * 4 + code1o;
78
79                code2o = amino_n[(int)s2[j][p2]];
80                code2u = amino_n[(int)s2[j][c2]];
81                if( code2o > 3 ) code2 = code1o = 36;
82                else if( code2u > 3 ) code2 = 36;
83                else code2 = 4 + code2u * 4 + code2o;
84
85
86//              fprintf( stderr, "%c%c-%c%c: %f\n", s1[i][p1], s1[i][c1], s2[j][p2], s2[j][c2], (float)ribosumdis[code1][code2] );
87
88                if( code1 == 36 || code2 == 36 )
89                        val += (float)n_dis[code1o][code2o] * eff1[i] * eff2[j];
90                else
91                        val += (float)ribosumdis[code1][code2] * eff1[i] * eff2[j];
92        }
93        return( val );
94}
95
96
97static void mccaskillextract( char **seq, char **nogap, int nseq, RNApair **pairprob, RNApair ***single, int **sgapmap, double *eff )
98{
99        int lgth;
100        int nogaplgth;
101        int i, j;
102        int left, right, adpos;
103        float prob;
104        static TLS int *pairnum;
105        RNApair *pt, *pt2;
106
107        lgth = strlen( seq[0] );
108        pairnum = calloc( lgth, sizeof( int ) );
109        for( i=0; i<lgth; i++ ) pairnum[i] = 0;
110
111        for( i=0; i<nseq; i++ )
112        {
113                nogaplgth = strlen( nogap[i] );
114                for( j=0; j<nogaplgth; j++ ) for( pt=single[i][j]; pt->bestpos!=-1; pt++ )
115                {
116                        left = sgapmap[i][j];
117                        right = sgapmap[i][pt->bestpos];
118                        prob = pt->bestscore;
119
120
121                        for( pt2=pairprob[left]; pt2->bestpos!=-1; pt2++ )
122                                if( pt2->bestpos == right ) break;
123
124//                      fprintf( stderr, "i,j=%d,%d, left=%d, right=%d, pt=%d, pt2->bestpos = %d\n", i, j, left, right, pt-single[i][j], pt2->bestpos );
125                        if( pt2->bestpos == -1 )
126                        {
127                                pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
128                                adpos = pairnum[left];
129                                pairnum[left]++;
130                                pairprob[left][adpos].bestscore = 0.0;
131                                pairprob[left][adpos].bestpos = right;
132                                pairprob[left][adpos+1].bestscore = -1.0;
133                                pairprob[left][adpos+1].bestpos = -1;
134                                pt2 = pairprob[left]+adpos;
135                        }
136                        else
137                                adpos = pt2-pairprob[left];
138
139                        pt2->bestscore += prob * eff[i];
140
141                        if( pt2->bestpos != right )
142                        {
143                                fprintf( stderr, "okashii!\n" );
144                                exit( 1 );
145                        }
146//                      fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
147//                      fprintf( stderr, "pairprob[0][0].bestpos=%d\n", pairprob[0][0].bestpos );
148//                      fprintf( stderr, "pairprob[0][0].bestscore=%f\n", pairprob[0][0].bestscore );
149                }
150        }
151
152//      fprintf( stderr, "before taikakuka\n" );
153        for( i=0; i<lgth; i++ ) for( j=0; j<pairnum[i]; j++ )
154        {
155                if( pairprob[i][j].bestpos > -1 )
156                {
157//                      pairprob[i][j].bestscore /= (float)nseq;
158//                      fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, pairprob[i][j].bestpos, pairprob[i][j].bestscore, seq[0][i], seq[0][pairprob[i][j].bestpos] );
159                }
160        }
161
162#if 0
163        for( i=0; i<lgth; i++ ) for( j=0; j<pairnum[i]; j++ )
164        {
165                right=pairprob[i][j].bestpos;
166                if( right < i ) continue;
167                fprintf( stderr, "no%d-%d, adding %d -> %d\n", i, j, right, i );
168                pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
169                pairprob[right][pairnum[right]].bestscore = pairprob[i][j].bestscore;
170                pairprob[right][pairnum[right]].bestpos = i;
171                pairnum[right]++;
172                pairprob[right][pairnum[right]].bestscore = -1.0;
173                pairprob[right][pairnum[right]].bestpos = -1;
174        }
175#endif
176
177        free( pairnum );
178
179}
180
181
182void rnaalifoldcall( char **seq, int nseq, RNApair **pairprob )
183{
184        int lgth;
185        int i;
186        static TLS int *order = NULL;
187        static TLS char **name = NULL;
188        char gett[1000];
189        FILE *fp;
190        int left, right, dumm;
191        float prob;
192        static TLS int pid;
193        static TLS char fnamein[100];
194        static TLS char cmd[1000];
195        static TLS int *pairnum;
196
197        lgth = strlen( seq[0] );
198        if( order == NULL )
199        {
200                pid = (int)getpid();
201                sprintf( fnamein, "/tmp/_rnaalifoldin.%d", pid );
202                order = AllocateIntVec( njob );
203                name = AllocateCharMtx( njob, 10 );
204                for( i=0; i<njob; i++ )
205                {
206                        order[i] = i;
207                        sprintf( name[i], "seq%d", i );
208                }
209        }
210        pairnum = calloc( lgth, sizeof( int ) );
211        for( i=0; i<lgth; i++ ) pairnum[i] = 0;
212
213        fp = fopen( fnamein, "w" );
214        if( !fp )
215        {
216                fprintf( stderr, "Cannot open /tmp/_rnaalifoldin\n" );
217                exit( 1 );
218        }
219        clustalout_pointer( fp, nseq, lgth, seq, name, NULL, NULL, order, 15 );
220        fclose( fp );
221
222        sprintf( cmd, "RNAalifold -p %s", fnamein );
223        system( cmd );
224
225        fp = fopen( "alifold.out", "r" );
226        if( !fp )
227        {
228                fprintf( stderr, "Cannot open /tmp/_rnaalifoldin\n" );
229                exit( 1 );
230        }
231
232#if 0
233        for( i=0; i<lgth; i++ ) // atode kesu
234        {
235                pairprob[i] = (RNApair *)realloc( pairprob[i], (2) * sizeof( RNApair ) ); // atode kesu
236                pairprob[i][1].bestscore = -1.0;
237                pairprob[i][1].bestpos = -1;
238        }
239#endif
240
241        while( 1 )
242        {
243                fgets( gett, 999, fp );
244                if( gett[0] == '(' ) break;
245                if( gett[0] == '{' ) break;
246                if( gett[0] == '.' ) break;
247                if( gett[0] == ',' ) break;
248                if( gett[0] != ' ' ) continue;
249
250                sscanf( gett, "%d %d %d %f", &left, &right, &dumm, &prob );
251                left--;
252                right--;
253
254
255#if 0
256                if( prob > 50.0 && prob > pairprob[left][0].bestscore )
257                {
258                        pairprob[left][0].bestscore = prob;
259                        pairprob[left][0].bestpos = right;
260#else
261                if( prob > 0.0 )
262                {
263                        pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
264                        pairprob[left][pairnum[left]].bestscore = prob / 100.0;
265                        pairprob[left][pairnum[left]].bestpos = right;
266                        pairnum[left]++;
267                        pairprob[left][pairnum[left]].bestscore = -1.0;
268                        pairprob[left][pairnum[left]].bestpos = -1;
269                        fprintf( stderr, "%d-%d, %f\n", left, right, prob );
270
271                        pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
272                        pairprob[right][pairnum[right]].bestscore = prob / 100.0;
273                        pairprob[right][pairnum[right]].bestpos = left;
274                        pairnum[right]++;
275                        pairprob[right][pairnum[right]].bestscore = -1.0;
276                        pairprob[right][pairnum[right]].bestpos = -1;
277                        fprintf( stderr, "%d-%d, %f\n", left, right, prob );
278#endif
279                }
280        }
281        fclose( fp );
282        sprintf( cmd, "rm -f %s", fnamein );
283        system( cmd ); 
284
285        for( i=0; i<lgth; i++ )
286        {
287                if( (right=pairprob[i][0].bestpos) > -1 )
288                {
289                        pairprob[right][0].bestpos = i;
290                        pairprob[right][0].bestscore = pairprob[i][0].bestscore;
291                }
292        }
293
294#if 0
295        for( i=0; i<lgth; i++ ) // atode kesu
296                if( pairprob[i][0].bestscore > -1 ) pairprob[i][0].bestscore = 1.0; // atode kesu
297#endif
298
299//      fprintf( stderr, "after taikakuka in rnaalifoldcall\n" );
300//      for( i=0; i<lgth; i++ )
301//      {
302//              fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, pairprob[i][0].bestpos, pairprob[i][0].bestscore, seq[0][i], seq[0][pairprob[i][0].bestpos] );
303//      }
304
305        free( pairnum );
306}
307
308
309static void utot( int n, int l, char **s )
310{
311        int i, j;
312        for( i=0; i<l; i++ )
313        {
314                for( j=0; j<n; j++ )
315                {
316                        if     ( s[j][i] == 'a' ) s[j][i] = 'a';
317                        else if( s[j][i] == 't' ) s[j][i] = 't';
318                        else if( s[j][i] == 'u' ) s[j][i] = 't';
319                        else if( s[j][i] == 'g' ) s[j][i] = 'g';
320                        else if( s[j][i] == 'c' ) s[j][i] = 'c';
321                        else if( s[j][i] == '-' ) s[j][i] = '-';
322                        else                                      s[j][i] = 'n';
323                }
324        }
325}
326
327
328void foldrna( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, float **impmtx, int *gapmap1, int *gapmap2, RNApair *additionalpair )
329{
330        int i, j;
331//      int ui, uj;
332//      int uiup, ujup;
333        int uido, ujdo;
334        static TLS char **useq1, **useq2;
335        static TLS char **oseq1, **oseq2, **oseq1r, **oseq2r, *odir1, *odir2;
336        static TLS RNApair **pairprob1, **pairprob2;
337        static TLS RNApair *pairpt1, *pairpt2;
338        int lgth1 = strlen( seq1[0] );
339        int lgth2 = strlen( seq2[0] );
340        static TLS float **impmtx2;
341        static TLS float **map;
342//      double lenfac;
343        float prob;
344        int **sgapmap1, **sgapmap2;
345//      char *nogapdum;
346        float **tbppmtx;
347
348
349//      fprintf( stderr, "nseq1=%d, lgth1=%d\n", nseq1, lgth1 );
350        useq1 = AllocateCharMtx( nseq1, lgth1+10 );
351        useq2 = AllocateCharMtx( nseq2, lgth2+10 );
352        oseq1 = AllocateCharMtx( nseq1, lgth1+10 );
353        oseq2 = AllocateCharMtx( nseq2, lgth2+10 );
354        oseq1r = AllocateCharMtx( nseq1, lgth1+10 );
355        oseq2r = AllocateCharMtx( nseq2, lgth2+10 );
356        odir1 = AllocateCharVec( lgth1+10 );
357        odir2 = AllocateCharVec( lgth2+10 );
358        sgapmap1 = AllocateIntMtx( nseq1, lgth1+1 );
359        sgapmap2 = AllocateIntMtx( nseq2, lgth2+1 );
360//      nogapdum = AllocateCharVec( MAX( lgth1, lgth2 ) );
361        pairprob1 = (RNApair **)calloc( lgth1, sizeof( RNApair *) );
362        pairprob2 = (RNApair **)calloc( lgth2, sizeof( RNApair *) );
363        map = AllocateFloatMtx( lgth1, lgth2 );
364        impmtx2 = AllocateFloatMtx( lgth1, lgth2 );
365        tbppmtx = AllocateFloatMtx( lgth1, lgth2 );
366
367        for( i=0; i<nseq1; i++ ) strcpy( useq1[i], seq1[i] );
368        for( i=0; i<nseq2; i++ ) strcpy( useq2[i], seq2[i] );
369        for( i=0; i<nseq1; i++ ) strcpy( oseq1[i], seq1[i] );
370        for( i=0; i<nseq2; i++ ) strcpy( oseq2[i], seq2[i] );
371
372        for( i=0; i<nseq1; i++ ) commongappick_record( 1, useq1+i, sgapmap1[i] );
373        for( i=0; i<nseq2; i++ ) commongappick_record( 1, useq2+i, sgapmap2[i] );
374
375        for( i=0; i<lgth1; i++ )
376        {
377                pairprob1[i] = (RNApair *)calloc( 1, sizeof( RNApair ) );
378                pairprob1[i][0].bestpos = -1;
379                pairprob1[i][0].bestscore = -1;
380        }
381        for( i=0; i<lgth2; i++ )
382        {
383                pairprob2[i] = (RNApair *)calloc( 1, sizeof( RNApair ) );
384                pairprob2[i][0].bestpos = -1;
385                pairprob2[i][0].bestscore = -1;
386        }
387
388        utot( nseq1, lgth1, oseq1 );
389        utot( nseq2, lgth2, oseq2 );
390
391//      fprintf( stderr, "folding group1\n" );
392//      rnalocal( oseq1, useq1, eff1, eff1, nseq1, nseq1, lgth1+10, pair1 );
393
394/* base-pairing probability of group 1 */
395        if( rnaprediction == 'r' )
396                rnaalifoldcall( oseq1, nseq1, pairprob1 );
397        else
398                mccaskillextract( oseq1, useq1, nseq1, pairprob1, grouprna1, sgapmap1, eff1 );
399
400
401//      fprintf( stderr, "folding group2\n" );
402//      rnalocal( oseq2, useq2, eff2, eff2, nseq2, nseq2, lgth2+10, pair2 );
403
404/* base-pairing probability of group 2 */
405        if( rnaprediction == 'r' )
406                rnaalifoldcall( oseq2, nseq2, pairprob2 );
407        else
408                mccaskillextract( oseq2, useq2, nseq2, pairprob2, grouprna2, sgapmap2, eff2 );
409
410
411
412#if 0
413        makerseq( oseq1, oseq1r, odir1, pairprob1, nseq1, lgth1 );
414        makerseq( oseq2, oseq2r, odir2, pairprob2, nseq2, lgth2 );
415
416        fprintf( stderr, "%s\n", odir2 );
417
418        for( i=0; i<nseq1; i++ )
419        {
420                fprintf( stdout, ">ori\n%s\n", oseq1[0] );
421                fprintf( stdout, ">rev\n%s\n", oseq1r[0] );
422        }
423#endif
424
425/* similarity score */
426        Lalignmm_hmout( oseq1, oseq2, eff1, eff2, nseq1, nseq2, 10000, NULL, NULL, NULL, NULL, map );
427
428        if( 1 )
429        {
430                if( RNAscoremtx == 'n' )
431                {
432                        for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
433                        {
434//                              impmtx2[i][j] = osoiaveragescore( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j ) * consweight_multi;
435                                impmtx2[i][j] = 0.0;
436                        }
437                }
438                else if( RNAscoremtx == 'r' )
439                {
440                        for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ ) 
441                        {
442                                tbppmtx[i][j] = 1.0;
443                                impmtx2[i][j] = 0.0;
444                        }
445                        for( i=0; i<lgth1; i++ ) for( pairpt1=pairprob1[i]; pairpt1->bestpos!=-1; pairpt1++ )
446                        {
447                                for( j=0; j<lgth2; j++ ) for( pairpt2=pairprob2[j]; pairpt2->bestpos!=-1; pairpt2++ )
448                                {
449                                        uido = pairpt1->bestpos;
450                                        ujdo = pairpt2->bestpos;
451                                        prob = pairpt1->bestscore * pairpt2->bestscore;
452                                        if( uido > -1  && ujdo > -1 )
453                                        {
454                                                if( uido > i && j > ujdo )
455                                                {
456                                                        impmtx2[i][j] += prob * pairedribosumscore53( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j, uido, ujdo ) * consweight_multi;
457                                                        tbppmtx[i][j] -= prob;
458                                                }
459                                                else if( i < uido && j < ujdo )
460                                                {
461                                                        impmtx2[i][j] += prob * pairedribosumscore35( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j, uido, ujdo ) * consweight_multi;
462                                                        tbppmtx[i][j] -= prob;
463                                                }
464                                        }
465                                }
466                        }
467       
468       
469                        for( i=0; i<lgth1; i++ )
470                        {
471                                for( j=0; j<lgth2; j++ )
472                                {
473                                        impmtx2[i][j] += tbppmtx[i][j] * singleribosumscore( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j ) * consweight_multi;
474                                }
475                        }
476                }
477
478
479/* four-way consistency */
480
481                for( i=0; i<lgth1; i++ ) for( pairpt1=pairprob1[i]; pairpt1->bestpos!=-1; pairpt1++ )
482                {
483
484//                      if( pairprob1[i] == NULL ) continue;
485
486                        for( j=0; j<lgth2; j++ ) for( pairpt2=pairprob2[j]; pairpt2->bestpos!=-1; pairpt2++ )
487                        {
488//                              fprintf( stderr, "i=%d, j=%d, pn1=%d, pn2=%d\n", i, j, pairpt1-pairprob1[i], pairpt2-pairprob2[j] );
489//                              if( pairprob2[j] == NULL ) continue;
490
491                                uido = pairpt1->bestpos;
492                                ujdo = pairpt2->bestpos;
493                                prob = pairpt1->bestscore * pairpt2->bestscore;
494//                              prob = 1.0;
495//                              fprintf( stderr, "i=%d->uido=%d, j=%d->ujdo=%d\n", i, uido, j, ujdo );
496
497//                              fprintf( stderr, "impmtx2[%d][%d] = %f\n", i, j, impmtx2[i][j] );
498
499//                              if( i < uido && j > ujdo ) continue;
500//                              if( i > uido && j < ujdo ) continue;
501
502
503//                              posdistj = abs( ujdo-j );
504
505//                              if( uido > -1 && ujdo > -1 ) 
506                                if( uido > -1 && ujdo > -1 && ( ( i > uido && j > ujdo ) || ( i < uido && j < ujdo ) ) )
507                                {
508                                        {
509                                                impmtx2[i][j] += MAX( 0, map[uido][ujdo] ) * consweight_rna * 600 * prob; // osoi
510                                        }
511                                }
512
513                        }
514                }
515                for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
516                {
517                        impmtx[i][j] += impmtx2[i][j];
518//                      fprintf( stderr, "fastathreshold=%f, consweight_multi=%f, consweight_rna=%f\n", fastathreshold, consweight_multi, consweight_rna );
519//                      impmtx[i][j] *= 0.5;
520                }
521
522//              impmtx[0][0] += 10000.0;
523//              impmtx[lgth1-1][lgth2-1] += 10000.0;
524
525
526
527#if 0
528                fprintf( stdout, "#impmtx2 = \n" );
529                for( i=0; i<lgth1; i++ )
530                {
531                        for( j=0; j<lgth2; j++ )
532                        {
533                                fprintf( stdout, "%d %d %f\n", i, j, impmtx2[i][j] );
534                        }
535                        fprintf( stdout, "\n" );
536                }
537                exit( 1 );
538#endif
539        }
540
541        FreeCharMtx( useq1 );
542        FreeCharMtx( useq2 );
543        FreeCharMtx( oseq1 );
544        FreeCharMtx( oseq2 );
545        FreeCharMtx( oseq1r );
546        FreeCharMtx( oseq2r );
547        free( odir1 );
548        free( odir2 );
549        FreeFloatMtx( impmtx2 );
550        FreeFloatMtx( map );
551        FreeIntMtx( sgapmap1 );
552        FreeIntMtx( sgapmap2 );
553        FreeFloatMtx( tbppmtx );
554
555        for( i=0; i<lgth1; i++ ) free( pairprob1[i] );
556        for( i=0; i<lgth2; i++ ) free( pairprob2[i] );
557        free( pairprob1 );
558        free( pairprob2 );
559}
560
Note: See TracBrowser for help on using the repository browser.