RNAlib-2.2.0-RC3
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 INLINE PRIVATE int backtrack_GQuad_IntLoop(int c,
147  int i,
148  int j,
149  int type,
150  short *S,
151  int *ggg,
152  int *index,
153  int *p,
154  int *q,
155  vrna_param_t *P);
156 
157 INLINE PRIVATE int backtrack_GQuad_IntLoop_L(int c,
158  int i,
159  int j,
160  int type,
161  short *S,
162  int **ggg,
163  int maxdist,
164  int *p,
165  int *q,
166  vrna_param_t *P);
167 
168 INLINE PRIVATE int
169 vrna_BT_gquad_mfe(vrna_fold_compound *vc,
170  int i,
171  int j,
172  bondT *bp_stack,
173  int *stack_count);
174 
175 INLINE PRIVATE int
176 vrna_BT_gquad_int(vrna_fold_compound *vc,
177  int i,
178  int j,
179  int en,
180  bondT *bp_stack,
181  int *stack_count);
182 
183 INLINE PRIVATE int
184 vrna_BT_gquad_mfe(vrna_fold_compound *vc,
185  int i,
186  int j,
187  bondT *bp_stack,
188  int *stack_count){
189 
190  /*
191  here we do some fancy stuff to backtrace the stacksize and linker lengths
192  of the g-quadruplex that should reside within position i,j
193  */
194  short *S;
195  int l[3], L, a;
196  vrna_param_t *P;
197 
198  P = vc->params;
199  S = vc->sequence_encoding2;
200  L = -1;
201 
202  get_gquad_pattern_mfe(S, i, j, P, &L, l);
203 
204  if(L != -1){
205  /* fill the G's of the quadruplex into base_pair2 */
206  for(a=0;a<L;a++){
207  bp_stack[++(*stack_count)].i = i+a;
208  bp_stack[(*stack_count)].j = i+a;
209  bp_stack[++(*stack_count)].i = i+L+l[0]+a;
210  bp_stack[(*stack_count)].j = i+L+l[0]+a;
211  bp_stack[++(*stack_count)].i = i+L+l[0]+L+l[1]+a;
212  bp_stack[(*stack_count)].j = i+L+l[0]+L+l[1]+a;
213  bp_stack[++(*stack_count)].i = i+L+l[0]+L+l[1]+L+l[2]+a;
214  bp_stack[(*stack_count)].j = i+L+l[0]+L+l[1]+L+l[2]+a;
215  }
216  return 1;
217  } else {
218  return 0;
219  }
220 }
221 
222 INLINE PRIVATE int
223 vrna_BT_gquad_int(vrna_fold_compound *vc,
224  int i,
225  int j,
226  int en,
227  bondT *bp_stack,
228  int *stack_count){
229 
230  int energy, dangles, *idx, ij, p, q, maxl, minl, c0, l1, *rtype, *ggg;
231  unsigned char type;
232  char *ptype;
233  short si, sj, *S, *S1;
234 
235  vrna_param_t *P;
236  vrna_md_t *md;
237 
238  idx = vc->jindx;
239  ij = idx[j] + i;
240  P = vc->params;
241  md = &(P->model_details);
242  ptype = vc->ptype;
243  rtype = &(md->rtype[0]);
244  type = rtype[(unsigned char)ptype[ij]];
245  S1 = vc->sequence_encoding;
246  S = vc->sequence_encoding2;
247  dangles = md->dangles;
248  si = S1[i + 1];
249  sj = S1[j - 1];
250  ggg = vc->matrices->ggg;
251  energy = 0;
252 
253  if(dangles == 2)
254  energy += P->mismatchI[type][si][sj];
255 
256  if(type > 2)
257  energy += P->TerminalAU;
258 
259  p = i + 1;
260  if(S1[p] == 3){
261  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
262  minl = j - i + p - MAXLOOP - 2;
263  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
264  minl = MAX2(c0, minl);
265  c0 = j - 3;
266  maxl = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
267  maxl = MIN2(c0, maxl);
268  for(q = minl; q < maxl; q++){
269  if(S[q] != 3) continue;
270  if(en == energy + ggg[idx[q] + p] + P->internal_loop[j - q - 1]){
271  return vrna_BT_gquad_mfe(vc, p, q, bp_stack, stack_count);
272  }
273  }
274  }
275  }
276 
277  for(p = i + 2;
278  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
279  p++){
280  l1 = p - i - 1;
281  if(l1>MAXLOOP) break;
282  if(S1[p] != 3) continue;
283  minl = j - i + p - MAXLOOP - 2;
284  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
285  minl = MAX2(c0, minl);
286  c0 = j - 1;
287  maxl = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
288  maxl = MIN2(c0, maxl);
289  for(q = minl; q < maxl; q++){
290  if(S1[q] != 3) continue;
291  if(en == energy + ggg[idx[q] + p] + P->internal_loop[l1 + j - q - 1]){
292  return vrna_BT_gquad_mfe(vc, p, q, bp_stack, stack_count);
293  }
294  }
295  }
296 
297  q = j - 1;
298  if(S1[q] == 3)
299  for(p = i + 4;
300  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
301  p++){
302  l1 = p - i - 1;
303  if(l1>MAXLOOP) break;
304  if(S1[p] != 3) continue;
305  if(en == energy + ggg[idx[q] + p] + P->internal_loop[l1]){
306  return vrna_BT_gquad_mfe(vc, p, q, bp_stack, stack_count);
307  }
308  }
309 
310  return 0;
311 }
312 
330 INLINE PRIVATE int backtrack_GQuad_IntLoop(int c,
331  int i,
332  int j,
333  int type,
334  short *S,
335  int *ggg,
336  int *index,
337  int *p,
338  int *q,
339  vrna_param_t *P){
340 
341  int energy, dangles, k, l, maxl, minl, c0, l1;
342  short si, sj;
343 
344  dangles = P->model_details.dangles;
345  si = S[i + 1];
346  sj = S[j - 1];
347  energy = 0;
348 
349  if(dangles == 2)
350  energy += P->mismatchI[type][si][sj];
351 
352  if(type > 2)
353  energy += P->TerminalAU;
354 
355  k = i + 1;
356  if(S[k] == 3){
357  if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
358  minl = j - i + k - MAXLOOP - 2;
359  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
360  minl = MAX2(c0, minl);
361  c0 = j - 3;
362  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
363  maxl = MIN2(c0, maxl);
364  for(l = minl; l < maxl; l++){
365  if(S[l] != 3) continue;
366  if(c == energy + ggg[index[l] + k] + P->internal_loop[j - l - 1]){
367  *p = k; *q = l;
368  return 1;
369  }
370  }
371  }
372  }
373 
374  for(k = i + 2;
375  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
376  k++){
377  l1 = k - i - 1;
378  if(l1>MAXLOOP) break;
379  if(S[k] != 3) continue;
380  minl = j - i + k - MAXLOOP - 2;
381  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
382  minl = MAX2(c0, minl);
383  c0 = j - 1;
384  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
385  maxl = MIN2(c0, maxl);
386  for(l = minl; l < maxl; l++){
387  if(S[l] != 3) continue;
388  if(c == energy + ggg[index[l] + k] + P->internal_loop[l1 + j - l - 1]){
389  *p = k; *q = l;
390  return 1;
391  }
392  }
393  }
394 
395  l = j - 1;
396  if(S[l] == 3)
397  for(k = i + 4;
398  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
399  k++){
400  l1 = k - i - 1;
401  if(l1>MAXLOOP) break;
402  if(S[k] != 3) continue;
403  if(c == energy + ggg[index[l] + k] + P->internal_loop[l1]){
404  *p = k; *q = l;
405  return 1;
406  }
407  }
408 
409  return 0;
410 }
411 
428 INLINE PRIVATE int backtrack_GQuad_IntLoop_L(int c,
429  int i,
430  int j,
431  int type,
432  short *S,
433  int **ggg,
434  int maxdist,
435  int *p,
436  int *q,
437  vrna_param_t *P){
438 
439  int energy, dangles, k, l, maxl, minl, c0, l1;
440  short si, sj;
441 
442  dangles = P->model_details.dangles;
443  si = S[i + 1];
444  sj = S[j - 1];
445  energy = 0;
446 
447  if(dangles == 2)
448  energy += P->mismatchI[type][si][sj];
449 
450  if(type > 2)
451  energy += P->TerminalAU;
452 
453  k = i + 1;
454  if(S[k] == 3){
455  if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
456  minl = j - i + k - MAXLOOP - 2;
457  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
458  minl = MAX2(c0, minl);
459  c0 = j - 3;
460  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
461  maxl = MIN2(c0, maxl);
462  for(l = minl; l < maxl; l++){
463  if(S[l] != 3) continue;
464  if(c == energy + ggg[k][l - k] + P->internal_loop[j - l - 1]){
465  *p = k; *q = l;
466  return 1;
467  }
468  }
469  }
470  }
471 
472  for(k = i + 2;
473  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
474  k++){
475  l1 = k - i - 1;
476  if(l1>MAXLOOP) break;
477  if(S[k] != 3) continue;
478  minl = j - i + k - MAXLOOP - 2;
479  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
480  minl = MAX2(c0, minl);
481  c0 = j - 1;
482  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
483  maxl = MIN2(c0, maxl);
484  for(l = minl; l < maxl; l++){
485  if(S[l] != 3) continue;
486  if(c == energy + ggg[k][l - k] + P->internal_loop[l1 + j - l - 1]){
487  *p = k; *q = l;
488  return 1;
489  }
490  }
491  }
492 
493  l = j - 1;
494  if(S[l] == 3)
495  for(k = i + 4;
496  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
497  k++){
498  l1 = k - i - 1;
499  if(l1>MAXLOOP) break;
500  if(S[k] != 3) continue;
501  if(c == energy + ggg[k][l - k] + P->internal_loop[l1]){
502  *p = k; *q = l;
503  return 1;
504  }
505  }
506 
507  return 0;
508 }
509 
510 INLINE PRIVATE
511 int
512 E_GQuad_IntLoop(int i,
513  int j,
514  int type,
515  short *S,
516  int *ggg,
517  int *index,
518  vrna_param_t *P){
519 
520  int energy, ge, dangles, p, q, l1, minq, maxq, c0;
521  short si, sj;
522 
523  dangles = P->model_details.dangles;
524  si = S[i + 1];
525  sj = S[j - 1];
526  energy = 0;
527 
528  if(dangles == 2)
529  energy += P->mismatchI[type][si][sj];
530 
531  if(type > 2)
532  energy += P->TerminalAU;
533 
534  ge = INF;
535 
536  p = i + 1;
537  if(S[p] == 3){
538  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
539  minq = j - i + p - MAXLOOP - 2;
540  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
541  minq = MAX2(c0, minq);
542  c0 = j - 3;
543  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
544  maxq = MIN2(c0, maxq);
545  for(q = minq; q < maxq; q++){
546  if(S[q] != 3) continue;
547  c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
548  ge = MIN2(ge, c0);
549  }
550  }
551  }
552 
553  for(p = i + 2;
554  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
555  p++){
556  l1 = p - i - 1;
557  if(l1>MAXLOOP) break;
558  if(S[p] != 3) continue;
559  minq = j - i + p - MAXLOOP - 2;
560  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
561  minq = MAX2(c0, minq);
562  c0 = j - 1;
563  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
564  maxq = MIN2(c0, maxq);
565  for(q = minq; q < maxq; q++){
566  if(S[q] != 3) continue;
567  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
568  ge = MIN2(ge, c0);
569  }
570  }
571 
572  q = j - 1;
573  if(S[q] == 3)
574  for(p = i + 4;
575  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
576  p++){
577  l1 = p - i - 1;
578  if(l1>MAXLOOP) break;
579  if(S[p] != 3) continue;
580  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
581  ge = MIN2(ge, c0);
582  }
583 
584 #if 0
585  /* here comes the additional stuff for the odd dangle models */
586  if(dangles % 1){
587  en1 = energy + P->dangle5[type][si];
588  en2 = energy + P->dangle5[type][sj];
589  en3 = energy + P->mismatchI[type][si][sj];
590 
591  /* first case with 5' dangle (i.e. j-1) onto enclosing pair */
592  p = i + 1;
593  if(S[p] == 3){
594  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
595  minq = j - i + p - MAXLOOP - 2;
596  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
597  minq = MAX2(c0, minq);
598  c0 = j - 4;
599  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
600  maxq = MIN2(c0, maxq);
601  for(q = minq; q < maxq; q++){
602  if(S[q] != 3) continue;
603  c0 = en1 + ggg[index[q] + p] + P->internal_loop[j - q - 1];
604  ge = MIN2(ge, c0);
605  }
606  }
607  }
608 
609  for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
610  l1 = p - i - 1;
611  if(l1>MAXLOOP) break;
612  if(S[p] != 3) continue;
613  minq = j - i + p - MAXLOOP - 2;
614  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
615  minq = MAX2(c0, minq);
616  c0 = j - 2;
617  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
618  maxq = MIN2(c0, maxq);
619  for(q = minq; q < maxq; q++){
620  if(S[q] != 3) continue;
621  c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
622  ge = MIN2(ge, c0);
623  }
624  }
625 
626  q = j - 2;
627  if(S[q] == 3)
628  for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
629  l1 = p - i - 1;
630  if(l1>MAXLOOP) break;
631  if(S[p] != 3) continue;
632  c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + 1];
633  ge = MIN2(ge, c0);
634  }
635 
636  /* second case with 3' dangle (i.e. i+1) onto enclosing pair */
637 
638  }
639 #endif
640  return ge;
641 }
642 
643 INLINE PRIVATE
644 int *
645 E_GQuad_IntLoop_exhaustive( int i,
646  int j,
647  int **p_p,
648  int **q_p,
649  int type,
650  short *S,
651  int *ggg,
652  int threshold,
653  int *index,
654  vrna_param_t *P){
655 
656  int energy, *ge, dangles, p, q, l1, minq, maxq, c0;
657  short si, sj;
658  int cnt = 0;
659 
660  dangles = P->model_details.dangles;
661  si = S[i + 1];
662  sj = S[j - 1];
663  energy = 0;
664 
665  if(dangles == 2)
666  energy += P->mismatchI[type][si][sj];
667 
668  if(type > 2)
669  energy += P->TerminalAU;
670 
671  /* guess how many gquads are possible in interval [i+1,j-1] */
672  *p_p = (int *)vrna_alloc(sizeof(int) * 256);
673  *q_p = (int *)vrna_alloc(sizeof(int) * 256);
674  ge = (int *)vrna_alloc(sizeof(int) * 256);
675 
676  p = i + 1;
677  if(S[p] == 3){
678  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
679  minq = j - i + p - MAXLOOP - 2;
680  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
681  minq = MAX2(c0, minq);
682  c0 = j - 3;
683  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
684  maxq = MIN2(c0, maxq);
685  for(q = minq; q < maxq; q++){
686  if(S[q] != 3) continue;
687  c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
688  if(c0 <= threshold){
689  ge[cnt] = energy + P->internal_loop[j - q - 1];
690  (*p_p)[cnt] = p;
691  (*q_p)[cnt++] = q;
692  }
693  }
694  }
695  }
696 
697  for(p = i + 2;
698  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
699  p++){
700  l1 = p - i - 1;
701  if(l1>MAXLOOP) break;
702  if(S[p] != 3) continue;
703  minq = j - i + p - MAXLOOP - 2;
704  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
705  minq = MAX2(c0, minq);
706  c0 = j - 1;
707  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
708  maxq = MIN2(c0, maxq);
709  for(q = minq; q < maxq; q++){
710  if(S[q] != 3) continue;
711  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
712  if(c0 <= threshold){
713  ge[cnt] = energy + P->internal_loop[l1 + j - q - 1];
714  (*p_p)[cnt] = p;
715  (*q_p)[cnt++] = q;
716  }
717  }
718  }
719 
720  q = j - 1;
721  if(S[q] == 3)
722  for(p = i + 4;
723  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
724  p++){
725  l1 = p - i - 1;
726  if(l1>MAXLOOP) break;
727  if(S[p] != 3) continue;
728  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
729  if(c0 <= threshold){
730  ge[cnt] = energy + P->internal_loop[l1];
731  (*p_p)[cnt] = p;
732  (*q_p)[cnt++] = q;
733  }
734  }
735 
736 
737  (*p_p)[cnt] = -1;
738 
739  return ge;
740 }
741 
742 INLINE PRIVATE
743 int
744 E_GQuad_IntLoop_L(int i,
745  int j,
746  int type,
747  short *S,
748  int **ggg,
749  int maxdist,
750  vrna_param_t *P){
751 
752  int energy, ge, dangles, p, q, l1, minq, maxq, c0;
753  short si, sj;
754 
755  dangles = P->model_details.dangles;
756  si = S[i + 1];
757  sj = S[j - 1];
758  energy = 0;
759 
760  if(dangles == 2)
761  energy += P->mismatchI[type][si][sj];
762 
763  if(type > 2)
764  energy += P->TerminalAU;
765 
766  ge = INF;
767 
768  p = i + 1;
769  if(S[p] == 3){
770  if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
771  minq = j - i + p - MAXLOOP - 2;
772  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
773  minq = MAX2(c0, minq);
774  c0 = j - 3;
775  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
776  maxq = MIN2(c0, maxq);
777  for(q = minq; q < maxq; q++){
778  if(S[q] != 3) continue;
779  c0 = energy + ggg[p][q-p] + P->internal_loop[j - q - 1];
780  ge = MIN2(ge, c0);
781  }
782  }
783  }
784 
785  for(p = i + 2;
786  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
787  p++){
788  l1 = p - i - 1;
789  if(l1>MAXLOOP) break;
790  if(S[p] != 3) continue;
791  minq = j - i + p - MAXLOOP - 2;
792  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
793  minq = MAX2(c0, minq);
794  c0 = j - 1;
795  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
796  maxq = MIN2(c0, maxq);
797  for(q = minq; q < maxq; q++){
798  if(S[q] != 3) continue;
799  c0 = energy + ggg[p][q - p] + P->internal_loop[l1 + j - q - 1];
800  ge = MIN2(ge, c0);
801  }
802  }
803 
804  q = j - 1;
805  if(S[q] == 3)
806  for(p = i + 4;
807  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
808  p++){
809  l1 = p - i - 1;
810  if(l1>MAXLOOP) break;
811  if(S[p] != 3) continue;
812  c0 = energy + ggg[p][q - p] + P->internal_loop[l1];
813  ge = MIN2(ge, c0);
814  }
815 
816  return ge;
817 }
818 
819 INLINE PRIVATE
820 FLT_OR_DBL
821 exp_E_GQuad_IntLoop(int i,
822  int j,
823  int type,
824  short *S,
825  FLT_OR_DBL *G,
826  int *index,
827  vrna_exp_param_t *pf){
828 
829  int k, l, minl, maxl, u, r;
830  FLT_OR_DBL q, qe, *expintern;
831  short si, sj;
832 
833  q = 0;
834  si = S[i + 1];
835  sj = S[j - 1];
836  qe = pf->expmismatchI[type][si][sj];
837  expintern = pf->expinternal;
838 
839  if(type > 2)
840  qe *= pf->expTermAU;
841 
842  k = i + 1;
843  if(S[k] == 3){
844  if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
845  minl = j - i + k - MAXLOOP - 2;
846  u = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
847  minl = MAX2(u, minl);
848  u = j - 3;
849  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
850  maxl = MIN2(u, maxl);
851  for(l = minl; l < maxl; l++){
852  if(S[l] != 3) continue;
853  if(G[index[k]-l] == 0.) continue;
854  q += qe * G[index[k]-l] * expintern[j - l - 1];
855  }
856  }
857  }
858 
859 
860  for(k = i + 2;
861  k <= j - VRNA_GQUAD_MIN_BOX_SIZE;
862  k++){
863  u = k - i - 1;
864  if(u > MAXLOOP) break;
865  if(S[k] != 3) continue;
866  minl = j - i + k - MAXLOOP - 2;
867  r = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
868  minl = MAX2(r, minl);
869  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
870  r = j - 1;
871  maxl = MIN2(r, maxl);
872  for(l = minl; l < maxl; l++){
873  if(S[l] != 3) continue;
874  if(G[index[k]-l] == 0.) continue;
875  q += qe * G[index[k]-l] * expintern[u + j - l - 1];
876  }
877  }
878 
879  l = j - 1;
880  if(S[l] == 3)
881  for(k = i + 4; k < j - VRNA_GQUAD_MIN_BOX_SIZE; k++){
882  u = k - i - 1;
883  if(u>MAXLOOP) break;
884  if(S[k] != 3) continue;
885  if(G[index[k]-l] == 0.) continue;
886  q += qe * G[index[k]-l] * expintern[u];
887  }
888 
889  return q;
890 }
891 
897 #endif
void * vrna_alloc(unsigned size)
Allocate space safely.
int * ggg
Energies of g-quadruplexes.
Definition: data_structures.h:432
#define MAXLOOP
Definition: energy_const.h:28
struct vrna_param_t * params
The precomputed free energy contributions for each type of loop.
Definition: data_structures.h:721
int * jindx
DP matrix accessor.
Definition: data_structures.h:725
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:428
The most basic data structure required by many functions throughout the RNAlib.
Definition: data_structures.h:698
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
vrna_mx_mfe_t * matrices
The MFE DP matrices.
Definition: data_structures.h:718
int dangles
Specifies the dangle model used in any energy evaluation (0,1,2 or 3)
Definition: model.h:172
char * ptype
Pair type array.
Definition: data_structures.h:749
Base pair.
Definition: data_structures.h:73
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
The data structure that contains the complete model details used throughout the calculations.
Definition: model.h:169
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:330
short * sequence_encoding
Numerical encoding of the input sequence.
Definition: data_structures.h:744
The datastructure that contains temperature scaled Boltzmann weights of the energy parameters...
Definition: params.h:86