RNAlib-2.1.9
gquad.h
Go to the documentation of this file.
1 #ifndef __VIENNA_RNA_PACKAGE_GQUAD_H__
2 #define __VIENNA_RNA_PACKAGE_GQUAD_H__
3 
4 #include "data_structures.h"
5 
6 #ifndef INLINE
7 #ifdef __GNUC__
8 # define INLINE inline
9 #else
10 # define INLINE
11 #endif
12 #endif
13 
20 int E_gquad(int L,
21  int l[3],
22  paramT *P);
23 
24 FLT_OR_DBL exp_E_gquad( int L,
25  int l[3],
26  pf_paramT *pf);
27 
28 int E_gquad_ali(int i,
29  int L,
30  int l[3],
31  const short **S,
32  int n_seq,
33  paramT *P);
34 
35 
36 void E_gquad_ali_en( int i,
37  int L,
38  int l[3],
39  const short **S,
40  int n_seq,
41  int en[2],
42  paramT *P);
43 
60 int *get_gquad_matrix(short *S, paramT *P);
61 
62 int *get_gquad_ali_matrix(short *S_cons,
63  short **S,
64  int n_seq,
65  paramT *P);
66 
67 FLT_OR_DBL *get_gquad_pf_matrix( short *S,
68  FLT_OR_DBL *scale,
69  pf_paramT *pf);
70 
71 int **get_gquad_L_matrix( short *S,
72  int start,
73  int maxdist,
74  int n,
75  int **g,
76  paramT *P);
77 
78 void get_gquad_pattern_mfe(short *S,
79  int i,
80  int j,
81  paramT *P,
82  int *L,
83  int l[3]);
84 
85 void
86 get_gquad_pattern_exhaustive( short *S,
87  int i,
88  int j,
89  paramT *P,
90  int *L,
91  int *l,
92  int threshold);
93 
94 void get_gquad_pattern_pf( short *S,
95  int i,
96  int j,
97  pf_paramT *pf,
98  int *L,
99  int l[3]);
100 
101 plist *get_plist_gquad_from_pr( short *S,
102  int gi,
103  int gj,
104  FLT_OR_DBL *G,
105  FLT_OR_DBL *probs,
106  FLT_OR_DBL *scale,
107  pf_paramT *pf);
108 plist *get_plist_gquad_from_pr_max(short *S,
109  int gi,
110  int gj,
111  FLT_OR_DBL *G,
112  FLT_OR_DBL *probs,
113  FLT_OR_DBL *scale,
114  int *L,
115  int l[3],
116  pf_paramT *pf);
117 
118 plist *get_plist_gquad_from_db( const char *structure,
119  float pr);
120 
121 int get_gquad_count(short *S,
122  int i,
123  int j);
124 
125 int get_gquad_layer_count(short *S,
126  int i,
127  int j);
128 
129 
140 int parse_gquad(const char *struc, int *L, int l[3]);
141 
142 
143 
161 INLINE PRIVATE int backtrack_GQuad_IntLoop(int c,
162  int i,
163  int j,
164  int type,
165  short *S,
166  int *ggg,
167  int *index,
168  int *p,
169  int *q,
170  paramT *P){
171 
172  int energy, dangles, k, l, maxl, minl, c0, l1;
173  short si, sj;
174 
175  dangles = P->model_details.dangles;
176  si = S[i + 1];
177  sj = S[j - 1];
178  energy = 0;
179 
180  if(dangles == 2)
181  energy += P->mismatchI[type][si][sj];
182 
183  if(type > 2)
184  energy += P->TerminalAU;
185 
186  k = i + 1;
187  if(S[k] == 3){
188  if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
189  minl = j - i + k - MAXLOOP - 2;
190  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
191  minl = MAX2(c0, minl);
192  c0 = j - 3;
193  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
194  maxl = MIN2(c0, maxl);
195  for(l = minl; l < maxl; l++){
196  if(S[l] != 3) continue;
197  if(c == energy + ggg[index[l] + k] + P->internal_loop[j - l - 1]){
198  *p = k; *q = l;
199  return 1;
200  }
201  }
202  }
203  }
204 
205  for(k = i + 2;
206  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
207  k++){
208  l1 = k - i - 1;
209  if(l1>MAXLOOP) break;
210  if(S[k] != 3) continue;
211  minl = j - i + k - MAXLOOP - 2;
212  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
213  minl = MAX2(c0, minl);
214  c0 = j - 1;
215  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
216  maxl = MIN2(c0, maxl);
217  for(l = minl; l < maxl; l++){
218  if(S[l] != 3) continue;
219  if(c == energy + ggg[index[l] + k] + P->internal_loop[l1 + j - l - 1]){
220  *p = k; *q = l;
221  return 1;
222  }
223  }
224  }
225 
226  l = j - 1;
227  if(S[l] == 3)
228  for(k = i + 4;
229  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
230  k++){
231  l1 = k - i - 1;
232  if(l1>MAXLOOP) break;
233  if(S[k] != 3) continue;
234  if(c == energy + ggg[index[l] + k] + P->internal_loop[l1]){
235  *p = k; *q = l;
236  return 1;
237  }
238  }
239 
240  return 0;
241 }
242 
259 INLINE PRIVATE int backtrack_GQuad_IntLoop_L(int c,
260  int i,
261  int j,
262  int type,
263  short *S,
264  int **ggg,
265  int maxdist,
266  int *p,
267  int *q,
268  paramT *P){
269 
270  int energy, dangles, k, l, maxl, minl, c0, l1;
271  short si, sj;
272 
273  dangles = P->model_details.dangles;
274  si = S[i + 1];
275  sj = S[j - 1];
276  energy = 0;
277 
278  if(dangles == 2)
279  energy += P->mismatchI[type][si][sj];
280 
281  if(type > 2)
282  energy += P->TerminalAU;
283 
284  k = i + 1;
285  if(S[k] == 3){
286  if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
287  minl = j - i + k - MAXLOOP - 2;
288  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
289  minl = MAX2(c0, minl);
290  c0 = j - 3;
291  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
292  maxl = MIN2(c0, maxl);
293  for(l = minl; l < maxl; l++){
294  if(S[l] != 3) continue;
295  if(c == energy + ggg[k][l - k] + P->internal_loop[j - l - 1]){
296  *p = k; *q = l;
297  return 1;
298  }
299  }
300  }
301  }
302 
303  for(k = i + 2;
304  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
305  k++){
306  l1 = k - i - 1;
307  if(l1>MAXLOOP) break;
308  if(S[k] != 3) continue;
309  minl = j - i + k - MAXLOOP - 2;
310  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
311  minl = MAX2(c0, minl);
312  c0 = j - 1;
313  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
314  maxl = MIN2(c0, maxl);
315  for(l = minl; l < maxl; l++){
316  if(S[l] != 3) continue;
317  if(c == energy + ggg[k][l - k] + P->internal_loop[l1 + j - l - 1]){
318  *p = k; *q = l;
319  return 1;
320  }
321  }
322  }
323 
324  l = j - 1;
325  if(S[l] == 3)
326  for(k = i + 4;
327  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
328  k++){
329  l1 = k - i - 1;
330  if(l1>MAXLOOP) break;
331  if(S[k] != 3) continue;
332  if(c == energy + ggg[k][l - k] + P->internal_loop[l1]){
333  *p = k; *q = l;
334  return 1;
335  }
336  }
337 
338  return 0;
339 }
340 
341 INLINE PRIVATE
342 int
343 E_GQuad_IntLoop(int i,
344  int j,
345  int type,
346  short *S,
347  int *ggg,
348  int *index,
349  paramT *P){
350 
351  int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
352  int c0, c1, c2, c3, up, d53, d5, d3;
353  short si, sj;
354 
355  dangles = P->model_details.dangles;
356  si = S[i + 1];
357  sj = S[j - 1];
358  energy = 0;
359 
360  if(dangles == 2)
361  energy += P->mismatchI[type][si][sj];
362 
363  if(type > 2)
364  energy += P->TerminalAU;
365 
366  ge = INF;
367 
368  p = i + 1;
369  if(S[p] == 3){
370  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
371  minq = j - i + p - MAXLOOP - 2;
372  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
373  minq = MAX2(c0, minq);
374  c0 = j - 3;
375  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
376  maxq = MIN2(c0, maxq);
377  for(q = minq; q < maxq; q++){
378  if(S[q] != 3) continue;
379  c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
380  ge = MIN2(ge, c0);
381  }
382  }
383  }
384 
385  for(p = i + 2;
386  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
387  p++){
388  l1 = p - i - 1;
389  if(l1>MAXLOOP) break;
390  if(S[p] != 3) continue;
391  minq = j - i + p - MAXLOOP - 2;
392  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
393  minq = MAX2(c0, minq);
394  c0 = j - 1;
395  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
396  maxq = MIN2(c0, maxq);
397  for(q = minq; q < maxq; q++){
398  if(S[q] != 3) continue;
399  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
400  ge = MIN2(ge, c0);
401  }
402  }
403 
404  q = j - 1;
405  if(S[q] == 3)
406  for(p = i + 4;
407  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
408  p++){
409  l1 = p - i - 1;
410  if(l1>MAXLOOP) break;
411  if(S[p] != 3) continue;
412  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
413  ge = MIN2(ge, c0);
414  }
415 
416 #if 0
417  /* here comes the additional stuff for the odd dangle models */
418  if(dangles % 1){
419  en1 = energy + P->dangle5[type][si];
420  en2 = energy + P->dangle5[type][sj];
421  en3 = energy + P->mismatchI[type][si][sj];
422 
423  /* first case with 5' dangle (i.e. j-1) onto enclosing pair */
424  p = i + 1;
425  if(S[p] == 3){
426  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
427  minq = j - i + p - MAXLOOP - 2;
428  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
429  minq = MAX2(c0, minq);
430  c0 = j - 4;
431  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
432  maxq = MIN2(c0, maxq);
433  for(q = minq; q < maxq; q++){
434  if(S[q] != 3) continue;
435  c0 = en1 + ggg[index[q] + p] + P->internal_loop[j - q - 1];
436  ge = MIN2(ge, c0);
437  }
438  }
439  }
440 
441  for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
442  l1 = p - i - 1;
443  if(l1>MAXLOOP) break;
444  if(S[p] != 3) continue;
445  minq = j - i + p - MAXLOOP - 2;
446  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
447  minq = MAX2(c0, minq);
448  c0 = j - 2;
449  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
450  maxq = MIN2(c0, maxq);
451  for(q = minq; q < maxq; q++){
452  if(S[q] != 3) continue;
453  c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
454  ge = MIN2(ge, c0);
455  }
456  }
457 
458  q = j - 2;
459  if(S[q] == 3)
460  for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
461  l1 = p - i - 1;
462  if(l1>MAXLOOP) break;
463  if(S[p] != 3) continue;
464  c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + 1];
465  ge = MIN2(ge, c0);
466  }
467 
468  /* second case with 3' dangle (i.e. i+1) onto enclosing pair */
469 
470  }
471 #endif
472  return ge;
473 }
474 
475 INLINE PRIVATE
476 int *
477 E_GQuad_IntLoop_exhaustive( int i,
478  int j,
479  int **p_p,
480  int **q_p,
481  int type,
482  short *S,
483  int *ggg,
484  int threshold,
485  int *index,
486  paramT *P){
487 
488  int energy, *ge, en1, en2, dangles, p, q, l1, minq, maxq;
489  int c0, c1, c2, c3, up, d53, d5, d3;
490  short si, sj;
491  int cnt = 0;
492 
493  dangles = P->model_details.dangles;
494  si = S[i + 1];
495  sj = S[j - 1];
496  energy = 0;
497 
498  if(dangles == 2)
499  energy += P->mismatchI[type][si][sj];
500 
501  if(type > 2)
502  energy += P->TerminalAU;
503 
504  /* guess how many gquads are possible in interval [i+1,j-1] */
505  *p_p = (int *)space(sizeof(int) * 256);
506  *q_p = (int *)space(sizeof(int) * 256);
507  ge = (int *)space(sizeof(int) * 256);
508 
509  p = i + 1;
510  if(S[p] == 3){
511  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
512  minq = j - i + p - MAXLOOP - 2;
513  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
514  minq = MAX2(c0, minq);
515  c0 = j - 3;
516  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
517  maxq = MIN2(c0, maxq);
518  for(q = minq; q < maxq; q++){
519  if(S[q] != 3) continue;
520  c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
521  if(c0 <= threshold){
522  ge[cnt] = energy + P->internal_loop[j - q - 1];
523  (*p_p)[cnt] = p;
524  (*q_p)[cnt++] = q;
525  }
526  }
527  }
528  }
529 
530  for(p = i + 2;
531  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
532  p++){
533  l1 = p - i - 1;
534  if(l1>MAXLOOP) break;
535  if(S[p] != 3) continue;
536  minq = j - i + p - MAXLOOP - 2;
537  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
538  minq = MAX2(c0, minq);
539  c0 = j - 1;
540  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
541  maxq = MIN2(c0, maxq);
542  for(q = minq; q < maxq; q++){
543  if(S[q] != 3) continue;
544  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
545  if(c0 <= threshold){
546  ge[cnt] = energy + P->internal_loop[l1 + j - q - 1];
547  (*p_p)[cnt] = p;
548  (*q_p)[cnt++] = q;
549  }
550  }
551  }
552 
553  q = j - 1;
554  if(S[q] == 3)
555  for(p = i + 4;
556  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
557  p++){
558  l1 = p - i - 1;
559  if(l1>MAXLOOP) break;
560  if(S[p] != 3) continue;
561  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
562  if(c0 <= threshold){
563  ge[cnt] = energy + P->internal_loop[l1];
564  (*p_p)[cnt] = p;
565  (*q_p)[cnt++] = q;
566  }
567  }
568 
569 
570  (*p_p)[cnt] = -1;
571 
572  return ge;
573 }
574 
575 INLINE PRIVATE
576 int
577 E_GQuad_IntLoop_L(int i,
578  int j,
579  int type,
580  short *S,
581  int **ggg,
582  int maxdist,
583  paramT *P){
584 
585  int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
586  int c0, c1, c2, c3, up, d53, d5, d3;
587  short si, sj;
588 
589  dangles = P->model_details.dangles;
590  si = S[i + 1];
591  sj = S[j - 1];
592  energy = 0;
593 
594  if(dangles == 2)
595  energy += P->mismatchI[type][si][sj];
596 
597  if(type > 2)
598  energy += P->TerminalAU;
599 
600  ge = INF;
601 
602  p = i + 1;
603  if(S[p] == 3){
604  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
605  minq = j - i + p - MAXLOOP - 2;
606  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
607  minq = MAX2(c0, minq);
608  c0 = j - 3;
609  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
610  maxq = MIN2(c0, maxq);
611  for(q = minq; q < maxq; q++){
612  if(S[q] != 3) continue;
613  c0 = energy + ggg[p][q-p] + P->internal_loop[j - q - 1];
614  ge = MIN2(ge, c0);
615  }
616  }
617  }
618 
619  for(p = i + 2;
620  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
621  p++){
622  l1 = p - i - 1;
623  if(l1>MAXLOOP) break;
624  if(S[p] != 3) continue;
625  minq = j - i + p - MAXLOOP - 2;
626  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
627  minq = MAX2(c0, minq);
628  c0 = j - 1;
629  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
630  maxq = MIN2(c0, maxq);
631  for(q = minq; q < maxq; q++){
632  if(S[q] != 3) continue;
633  c0 = energy + ggg[p][q - p] + P->internal_loop[l1 + j - q - 1];
634  ge = MIN2(ge, c0);
635  }
636  }
637 
638  q = j - 1;
639  if(S[q] == 3)
640  for(p = i + 4;
641  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
642  p++){
643  l1 = p - i - 1;
644  if(l1>MAXLOOP) break;
645  if(S[p] != 3) continue;
646  c0 = energy + ggg[p][q - p] + P->internal_loop[l1];
647  ge = MIN2(ge, c0);
648  }
649 
650  return ge;
651 }
652 
653 INLINE PRIVATE
654 FLT_OR_DBL
655 exp_E_GQuad_IntLoop(int i,
656  int j,
657  int type,
658  short *S,
659  FLT_OR_DBL *G,
660  int *index,
661  pf_paramT *pf){
662 
663  int k, l, minl, maxl, u, r;
664  FLT_OR_DBL q, qe, *expintern;
665  short si, sj;
666 
667  q = 0;
668  si = S[i + 1];
669  sj = S[j - 1];
670  qe = pf->expmismatchI[type][si][sj];
671  expintern = pf->expinternal;
672 
673  if(type > 2)
674  qe *= pf->expTermAU;
675 
676  k = i + 1;
677  if(S[k] == 3){
678  if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
679  minl = j - i + k - MAXLOOP - 2;
680  u = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
681  minl = MAX2(u, minl);
682  u = j - 3;
683  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
684  maxl = MIN2(u, maxl);
685  for(l = minl; l < maxl; l++){
686  if(S[l] != 3) continue;
687  if(G[index[k]-l] == 0.) continue;
688  q += qe * G[index[k]-l] * expintern[j - l - 1];
689  }
690  }
691  }
692 
693 
694  for(k = i + 2;
695  k <= j - VRNA_GQUAD_MIN_BOX_SIZE;
696  k++){
697  u = k - i - 1;
698  if(u > MAXLOOP) break;
699  if(S[k] != 3) continue;
700  minl = j - i + k - MAXLOOP - 2;
701  r = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
702  minl = MAX2(r, minl);
703  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
704  r = j - 1;
705  maxl = MIN2(r, maxl);
706  for(l = minl; l < maxl; l++){
707  if(S[l] != 3) continue;
708  if(G[index[k]-l] == 0.) continue;
709  q += qe * G[index[k]-l] * expintern[u + j - l - 1];
710  }
711  }
712 
713  l = j - 1;
714  if(S[l] == 3)
715  for(k = i + 4; k < j - VRNA_GQUAD_MIN_BOX_SIZE; k++){
716  u = k - i - 1;
717  if(u>MAXLOOP) break;
718  if(S[k] != 3) continue;
719  if(G[index[k]-l] == 0.) continue;
720  q += qe * G[index[k]-l] * expintern[u];
721  }
722 
723  return q;
724 }
725 
726 #endif