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