RNAlib-2.4.14
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 <ViennaRNA/datastructures/basic.h>
5 #include <ViennaRNA/fold_compound.h>
6 #include <ViennaRNA/params/basic.h>
7 
8 #ifndef INLINE
9 #ifdef __GNUC__
10 # define INLINE inline
11 #else
12 # define INLINE
13 #endif
14 #endif
15 
30 int E_gquad(int L,
31  int l[3],
32  vrna_param_t *P);
33 
34 
35 FLT_OR_DBL exp_E_gquad(int L,
36  int l[3],
37  vrna_exp_param_t *pf);
38 
39 
40 void E_gquad_ali_en(int i,
41  int L,
42  int l[3],
43  const short **S,
44  unsigned int **a2s,
45  unsigned int n_seq,
46  vrna_param_t *P,
47  int en[2]);
48 
49 
66 int *get_gquad_matrix(short *S,
67  vrna_param_t *P);
68 
69 
70 int *get_gquad_ali_matrix(short *S_cons,
71  short **S,
72  unsigned int **a2s,
73  int n_seq,
74  vrna_param_t *P);
75 
76 
77 FLT_OR_DBL *get_gquad_pf_matrix(short *S,
78  FLT_OR_DBL *scale,
79  vrna_exp_param_t *pf);
80 
81 
82 FLT_OR_DBL *get_gquad_pf_matrix_comparative(short *S_cons,
83  short **S,
84  unsigned int **a2s,
85  FLT_OR_DBL *scale,
86  unsigned int n_seq,
87  vrna_exp_param_t *pf);
88 
89 
90 int **get_gquad_L_matrix(short *S,
91  int start,
92  int maxdist,
93  int n,
94  int **g,
95  vrna_param_t *P);
96 
97 
98 void vrna_gquad_mx_local_update(vrna_fold_compound_t *vc,
99  int start);
100 
101 
102 void get_gquad_pattern_mfe(short *S,
103  int i,
104  int j,
105  vrna_param_t *P,
106  int *L,
107  int l[3]);
108 
109 
110 void
111 get_gquad_pattern_exhaustive(short *S,
112  int i,
113  int j,
114  vrna_param_t *P,
115  int *L,
116  int *l,
117  int threshold);
118 
119 
120 void get_gquad_pattern_pf(short *S,
121  int i,
122  int j,
123  vrna_exp_param_t *pf,
124  int *L,
125  int l[3]);
126 
127 
128 plist *get_plist_gquad_from_pr(short *S,
129  int gi,
130  int gj,
131  FLT_OR_DBL *G,
132  FLT_OR_DBL *probs,
133  FLT_OR_DBL *scale,
134  vrna_exp_param_t *pf);
135 
136 
137 plist *get_plist_gquad_from_pr_max(short *S,
138  int gi,
139  int gj,
140  FLT_OR_DBL *G,
141  FLT_OR_DBL *probs,
142  FLT_OR_DBL *scale,
143  int *L,
144  int l[3],
145  vrna_exp_param_t *pf);
146 
147 
148 plist *get_plist_gquad_from_db(const char *structure,
149  float pr);
150 
151 
152 plist *
153 vrna_get_plist_gquad_from_pr(vrna_fold_compound_t *fc,
154  int gi,
155  int gj);
156 
157 
158 plist *
159 vrna_get_plist_gquad_from_pr_max(vrna_fold_compound_t *fc,
160  int gi,
161  int gj,
162  int *Lmax,
163  int lmax[3]);
164 
165 
166 int get_gquad_count(short *S,
167  int i,
168  int j);
169 
170 
171 int get_gquad_layer_count(short *S,
172  int i,
173  int j);
174 
175 
176 void get_gquad_pattern_mfe_ali(short **S,
177  unsigned int **a2s,
178  short *S_cons,
179  int n_seq,
180  int i,
181  int j,
182  vrna_param_t *P,
183  int *L,
184  int l[3]);
185 
186 
197 int parse_gquad(const char *struc,
198  int *L,
199  int l[3]);
200 
201 
202 INLINE PRIVATE int backtrack_GQuad_IntLoop(int c,
203  int i,
204  int j,
205  int type,
206  short *S,
207  int *ggg,
208  int *index,
209  int *p,
210  int *q,
211  vrna_param_t *P);
212 
213 
214 INLINE PRIVATE int backtrack_GQuad_IntLoop_comparative(int c,
215  int i,
216  int j,
217  unsigned int *type,
218  short *S_cons,
219  short **S5,
220  short **S3,
221  unsigned int **a2s,
222  int *ggg,
223  int *index,
224  int *p,
225  int *q,
226  int n_seq,
227  vrna_param_t *P);
228 
229 
230 INLINE PRIVATE int backtrack_GQuad_IntLoop_L(int c,
231  int i,
232  int j,
233  int type,
234  short *S,
235  int **ggg,
236  int maxdist,
237  int *p,
238  int *q,
239  vrna_param_t *P);
240 
241 
242 PRIVATE INLINE int
243 vrna_BT_gquad_int(vrna_fold_compound_t *vc,
244  int i,
245  int j,
246  int en,
247  vrna_bp_stack_t *bp_stack,
248  int *stack_count);
249 
250 
251 PRIVATE INLINE int
252 vrna_BT_gquad_mfe(vrna_fold_compound_t *vc,
253  int i,
254  int j,
255  vrna_bp_stack_t *bp_stack,
256  int *stack_count)
257 {
258  /*
259  * here we do some fancy stuff to backtrace the stacksize and linker lengths
260  * of the g-quadruplex that should reside within position i,j
261  */
262  short *S;
263  int l[3], L, a, n_seq;
264  vrna_param_t *P;
265 
266  if (vc) {
267  P = vc->params;
268  switch (vc->type) {
269  case VRNA_FC_TYPE_SINGLE:
270  S = vc->sequence_encoding2;
271  L = -1;
272 
273  get_gquad_pattern_mfe(S, i, j, P, &L, l);
274  break;
275 
276  case VRNA_FC_TYPE_COMPARATIVE:
277  n_seq = vc->n_seq;
278  L = -1;
279  get_gquad_pattern_mfe_ali(vc->S, vc->a2s, vc->S_cons, n_seq, i, j, P, &L, l);
280  break;
281  }
282 
283  if (L != -1) {
284  /* fill the G's of the quadruplex into base_pair2 */
285  for (a = 0; a < L; a++) {
286  bp_stack[++(*stack_count)].i = i + a;
287  bp_stack[(*stack_count)].j = i + a;
288  bp_stack[++(*stack_count)].i = i + L + l[0] + a;
289  bp_stack[(*stack_count)].j = i + L + l[0] + a;
290  bp_stack[++(*stack_count)].i = i + L + l[0] + L + l[1] + a;
291  bp_stack[(*stack_count)].j = i + L + l[0] + L + l[1] + a;
292  bp_stack[++(*stack_count)].i = i + L + l[0] + L + l[1] + L + l[2] + a;
293  bp_stack[(*stack_count)].j = i + L + l[0] + L + l[1] + L + l[2] + a;
294  }
295  return 1;
296  } else {
297  return 0;
298  }
299  }
300 
301  return 0;
302 }
303 
304 
305 PRIVATE INLINE int
306 vrna_BT_gquad_int(vrna_fold_compound_t *vc,
307  int i,
308  int j,
309  int en,
310  vrna_bp_stack_t *bp_stack,
311  int *stack_count)
312 {
313  int energy, dangles, *idx, ij, p, q, maxl, minl, c0, l1, *ggg;
314  unsigned char type;
315  char *ptype;
316  short si, sj, *S, *S1;
317 
318  vrna_param_t *P;
319  vrna_md_t *md;
320 
321  idx = vc->jindx;
322  ij = idx[j] + i;
323  P = vc->params;
324  md = &(P->model_details);
325  ptype = vc->ptype;
326  type = (unsigned char)ptype[ij];
327  S1 = vc->sequence_encoding;
328  S = vc->sequence_encoding2;
329  dangles = md->dangles;
330  si = S1[i + 1];
331  sj = S1[j - 1];
332  ggg = vc->matrices->ggg;
333  energy = 0;
334 
335  if (dangles == 2)
336  energy += P->mismatchI[type][si][sj];
337 
338  if (type > 2)
339  energy += P->TerminalAU;
340 
341  p = i + 1;
342  if (S1[p] == 3) {
343  if (p < j - VRNA_GQUAD_MIN_BOX_SIZE) {
344  minl = j - i + p - MAXLOOP - 2;
345  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
346  minl = MAX2(c0, minl);
347  c0 = j - 3;
348  maxl = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
349  maxl = MIN2(c0, maxl);
350  for (q = minl; q < maxl; q++) {
351  if (S[q] != 3)
352  continue;
353 
354  if (en == energy + ggg[idx[q] + p] + P->internal_loop[j - q - 1])
355  return vrna_BT_gquad_mfe(vc, p, q, bp_stack, stack_count);
356  }
357  }
358  }
359 
360  for (p = i + 2;
361  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
362  p++) {
363  l1 = p - i - 1;
364  if (l1 > MAXLOOP)
365  break;
366 
367  if (S1[p] != 3)
368  continue;
369 
370  minl = j - i + p - MAXLOOP - 2;
371  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
372  minl = MAX2(c0, minl);
373  c0 = j - 1;
374  maxl = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
375  maxl = MIN2(c0, maxl);
376  for (q = minl; q < maxl; q++) {
377  if (S1[q] != 3)
378  continue;
379 
380  if (en == energy + ggg[idx[q] + p] + P->internal_loop[l1 + j - q - 1])
381  return vrna_BT_gquad_mfe(vc, p, q, bp_stack, stack_count);
382  }
383  }
384 
385  q = j - 1;
386  if (S1[q] == 3)
387  for (p = i + 4;
388  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
389  p++) {
390  l1 = p - i - 1;
391  if (l1 > MAXLOOP)
392  break;
393 
394  if (S1[p] != 3)
395  continue;
396 
397  if (en == energy + ggg[idx[q] + p] + P->internal_loop[l1])
398  return vrna_BT_gquad_mfe(vc, p, q, bp_stack, stack_count);
399  }
400 
401  return 0;
402 }
403 
404 
422 INLINE PRIVATE int
423 backtrack_GQuad_IntLoop(int c,
424  int i,
425  int j,
426  int type,
427  short *S,
428  int *ggg,
429  int *index,
430  int *p,
431  int *q,
432  vrna_param_t *P)
433 {
434  int energy, dangles, k, l, maxl, minl, c0, l1;
435  short si, sj;
436 
437  dangles = P->model_details.dangles;
438  si = S[i + 1];
439  sj = S[j - 1];
440  energy = 0;
441 
442  if (dangles == 2)
443  energy += P->mismatchI[type][si][sj];
444 
445  if (type > 2)
446  energy += P->TerminalAU;
447 
448  k = i + 1;
449  if (S[k] == 3) {
450  if (k < j - VRNA_GQUAD_MIN_BOX_SIZE) {
451  minl = j - i + k - MAXLOOP - 2;
452  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
453  minl = MAX2(c0, minl);
454  c0 = j - 3;
455  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
456  maxl = MIN2(c0, maxl);
457  for (l = minl; l < maxl; l++) {
458  if (S[l] != 3)
459  continue;
460 
461  if (c == energy + ggg[index[l] + k] + P->internal_loop[j - l - 1]) {
462  *p = k;
463  *q = l;
464  return 1;
465  }
466  }
467  }
468  }
469 
470  for (k = i + 2;
471  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
472  k++) {
473  l1 = k - i - 1;
474  if (l1 > MAXLOOP)
475  break;
476 
477  if (S[k] != 3)
478  continue;
479 
480  minl = j - i + k - MAXLOOP - 2;
481  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
482  minl = MAX2(c0, minl);
483  c0 = j - 1;
484  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
485  maxl = MIN2(c0, maxl);
486  for (l = minl; l < maxl; l++) {
487  if (S[l] != 3)
488  continue;
489 
490  if (c == energy + ggg[index[l] + k] + P->internal_loop[l1 + j - l - 1]) {
491  *p = k;
492  *q = l;
493  return 1;
494  }
495  }
496  }
497 
498  l = j - 1;
499  if (S[l] == 3)
500  for (k = i + 4;
501  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
502  k++) {
503  l1 = k - i - 1;
504  if (l1 > MAXLOOP)
505  break;
506 
507  if (S[k] != 3)
508  continue;
509 
510  if (c == energy + ggg[index[l] + k] + P->internal_loop[l1]) {
511  *p = k;
512  *q = l;
513  return 1;
514  }
515  }
516 
517  return 0;
518 }
519 
520 
521 INLINE PRIVATE int
522 backtrack_GQuad_IntLoop_comparative(int c,
523  int i,
524  int j,
525  unsigned int *type,
526  short *S_cons,
527  short **S5,
528  short **S3,
529  unsigned int **a2s,
530  int *ggg,
531  int *index,
532  int *p,
533  int *q,
534  int n_seq,
535  vrna_param_t *P)
536 {
537  int energy, dangles, k, l, maxl, minl, c0, l1, ss, tt, u1, u2, eee;
538 
539  dangles = P->model_details.dangles;
540  energy = 0;
541 
542  for (ss = 0; ss < n_seq; ss++) {
543  tt = type[ss];
544  if (tt == 0)
545  tt = 7;
546 
547  if (dangles == 2)
548  energy += P->mismatchI[tt][S3[ss][i]][S5[ss][j]];
549 
550  if (tt > 2)
551  energy += P->TerminalAU;
552  }
553 
554  k = i + 1;
555  if (S_cons[k] == 3) {
556  if (k < j - VRNA_GQUAD_MIN_BOX_SIZE) {
557  minl = j - i + k - MAXLOOP - 2;
558  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
559  minl = MAX2(c0, minl);
560  c0 = j - 3;
561  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
562  maxl = MIN2(c0, maxl);
563  for (l = minl; l < maxl; l++) {
564  if (S_cons[l] != 3)
565  continue;
566 
567  eee = 0;
568 
569  for (ss = 0; ss < n_seq; ss++) {
570  u1 = a2s[ss][j - 1] - a2s[ss][l];
571  eee += P->internal_loop[u1];
572  }
573 
574  if (c == energy + ggg[index[l] + k] + eee) {
575  *p = k;
576  *q = l;
577  return 1;
578  }
579  }
580  }
581  }
582 
583  for (k = i + 2;
584  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
585  k++) {
586  l1 = k - i - 1;
587  if (l1 > MAXLOOP)
588  break;
589 
590  if (S_cons[k] != 3)
591  continue;
592 
593  minl = j - i + k - MAXLOOP - 2;
594  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
595  minl = MAX2(c0, minl);
596  c0 = j - 1;
597  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
598  maxl = MIN2(c0, maxl);
599  for (l = minl; l < maxl; l++) {
600  if (S_cons[l] != 3)
601  continue;
602 
603  eee = 0;
604 
605  for (ss = 0; ss < n_seq; ss++) {
606  u1 = a2s[ss][k - 1] - a2s[ss][i];
607  u2 = a2s[ss][j - 1] - a2s[ss][l];
608  eee += P->internal_loop[u1 + u2];
609  }
610 
611  if (c == energy + ggg[index[l] + k] + eee) {
612  *p = k;
613  *q = l;
614  return 1;
615  }
616  }
617  }
618 
619  l = j - 1;
620  if (S_cons[l] == 3)
621  for (k = i + 4;
622  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
623  k++) {
624  l1 = k - i - 1;
625  if (l1 > MAXLOOP)
626  break;
627 
628  if (S_cons[k] != 3)
629  continue;
630 
631  eee = 0;
632 
633  for (ss = 0; ss < n_seq; ss++) {
634  u1 = a2s[ss][k - 1] - a2s[ss][i];
635  eee += P->internal_loop[u1];
636  }
637 
638  if (c == energy + ggg[index[l] + k] + eee) {
639  *p = k;
640  *q = l;
641  return 1;
642  }
643  }
644 
645  return 0;
646 }
647 
648 
665 INLINE PRIVATE int
666 backtrack_GQuad_IntLoop_L(int c,
667  int i,
668  int j,
669  int type,
670  short *S,
671  int **ggg,
672  int maxdist,
673  int *p,
674  int *q,
675  vrna_param_t *P)
676 {
677  int energy, dangles, k, l, maxl, minl, c0, l1;
678  short si, sj;
679 
680  dangles = P->model_details.dangles;
681  si = S[i + 1];
682  sj = S[j - 1];
683  energy = 0;
684 
685  if (dangles == 2)
686  energy += P->mismatchI[type][si][sj];
687 
688  if (type > 2)
689  energy += P->TerminalAU;
690 
691  k = i + 1;
692  if (S[k] == 3) {
693  if (k < j - VRNA_GQUAD_MIN_BOX_SIZE) {
694  minl = j - i + k - MAXLOOP - 2;
695  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
696  minl = MAX2(c0, minl);
697  c0 = j - 3;
698  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
699  maxl = MIN2(c0, maxl);
700  for (l = minl; l < maxl; l++) {
701  if (S[l] != 3)
702  continue;
703 
704  if (c == energy + ggg[k][l - k] + P->internal_loop[j - l - 1]) {
705  *p = k;
706  *q = l;
707  return 1;
708  }
709  }
710  }
711  }
712 
713  for (k = i + 2;
714  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
715  k++) {
716  l1 = k - i - 1;
717  if (l1 > MAXLOOP)
718  break;
719 
720  if (S[k] != 3)
721  continue;
722 
723  minl = j - i + k - MAXLOOP - 2;
724  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
725  minl = MAX2(c0, minl);
726  c0 = j - 1;
727  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
728  maxl = MIN2(c0, maxl);
729  for (l = minl; l < maxl; l++) {
730  if (S[l] != 3)
731  continue;
732 
733  if (c == energy + ggg[k][l - k] + P->internal_loop[l1 + j - l - 1]) {
734  *p = k;
735  *q = l;
736  return 1;
737  }
738  }
739  }
740 
741  l = j - 1;
742  if (S[l] == 3)
743  for (k = i + 4;
744  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
745  k++) {
746  l1 = k - i - 1;
747  if (l1 > MAXLOOP)
748  break;
749 
750  if (S[k] != 3)
751  continue;
752 
753  if (c == energy + ggg[k][l - k] + P->internal_loop[l1]) {
754  *p = k;
755  *q = l;
756  return 1;
757  }
758  }
759 
760  return 0;
761 }
762 
763 
764 INLINE PRIVATE int
765 backtrack_GQuad_IntLoop_L_comparative(int c,
766  int i,
767  int j,
768  unsigned int *type,
769  short *S_cons,
770  short **S5,
771  short **S3,
772  unsigned int **a2s,
773  int **ggg,
774  int *p,
775  int *q,
776  int n_seq,
777  vrna_param_t *P)
778 {
779  /*
780  * The case that is handled here actually resembles something like
781  * an interior loop where the enclosing base pair is of regular
782  * kind and the enclosed pair is not a canonical one but a g-quadruplex
783  * that should then be decomposed further...
784  */
785  int mm, dangle_model, k, l, maxl, minl, c0, l1, ss, tt, eee, u1, u2;
786 
787  dangle_model = P->model_details.dangles;
788 
789  mm = 0;
790  for (ss = 0; ss < n_seq; ss++) {
791  tt = type[ss];
792 
793  if (dangle_model == 2)
794  mm += P->mismatchI[tt][S3[ss][i]][S5[ss][j]];
795 
796  if (tt > 2)
797  mm += P->TerminalAU;
798  }
799 
800  for (k = i + 2;
801  k < j - VRNA_GQUAD_MIN_BOX_SIZE;
802  k++) {
803  if (S_cons[k] != 3)
804  continue;
805 
806  l1 = k - i - 1;
807  if (l1 > MAXLOOP)
808  break;
809 
810  minl = j - i + k - MAXLOOP - 2;
811  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
812  minl = MAX2(c0, minl);
813  c0 = j - 1;
814  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
815  maxl = MIN2(c0, maxl);
816  for (l = minl; l < maxl; l++) {
817  if (S_cons[l] != 3)
818  continue;
819 
820  eee = 0;
821 
822  for (ss = 0; ss < n_seq; ss++) {
823  u1 = a2s[ss][k - 1] - a2s[ss][i];
824  u2 = a2s[ss][j - 1] - a2s[ss][l];
825  eee += P->internal_loop[u1 + u2];
826  }
827 
828  c0 = mm +
829  ggg[k][l - k] +
830  eee;
831 
832  if (c == c0) {
833  *p = k;
834  *q = l;
835  return 1;
836  }
837  }
838  }
839  k = i + 1;
840  if (S_cons[k] == 3) {
841  if (k < j - VRNA_GQUAD_MIN_BOX_SIZE) {
842  minl = j - i + k - MAXLOOP - 2;
843  c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
844  minl = MAX2(c0, minl);
845  c0 = j - 3;
846  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
847  maxl = MIN2(c0, maxl);
848  for (l = minl; l < maxl; l++) {
849  if (S_cons[l] != 3)
850  continue;
851 
852  eee = 0;
853 
854  for (ss = 0; ss < n_seq; ss++) {
855  u1 = a2s[ss][j - 1] - a2s[ss][l];
856  eee += P->internal_loop[u1];
857  }
858 
859  if (c == mm + ggg[k][l - k] + eee) {
860  *p = k;
861  *q = l;
862  return 1;
863  }
864  }
865  }
866  }
867 
868  l = j - 1;
869  if (S_cons[l] == 3) {
870  for (k = i + 4; k < j - VRNA_GQUAD_MIN_BOX_SIZE; k++) {
871  l1 = k - i - 1;
872  if (l1 > MAXLOOP)
873  break;
874 
875  if (S_cons[k] != 3)
876  continue;
877 
878  eee = 0;
879 
880  for (ss = 0; ss < n_seq; ss++) {
881  u1 = a2s[ss][k - 1] - a2s[ss][i];
882  eee += P->internal_loop[u1];
883  }
884 
885  if (c == mm + ggg[k][l - k] + eee) {
886  *p = k;
887  *q = l;
888  return 1;
889  }
890  }
891  }
892 
893  return 0;
894 }
895 
896 
897 PRIVATE INLINE
898 int
899 E_GQuad_IntLoop(int i,
900  int j,
901  int type,
902  short *S,
903  int *ggg,
904  int *index,
905  vrna_param_t *P)
906 {
907  int energy, ge, dangles, p, q, l1, minq, maxq, c0;
908  short si, sj;
909 
910  dangles = P->model_details.dangles;
911  si = S[i + 1];
912  sj = S[j - 1];
913  energy = 0;
914 
915  if (dangles == 2)
916  energy += P->mismatchI[type][si][sj];
917 
918  if (type > 2)
919  energy += P->TerminalAU;
920 
921  ge = INF;
922 
923  p = i + 1;
924  if (S[p] == 3) {
925  if (p < j - VRNA_GQUAD_MIN_BOX_SIZE) {
926  minq = j - i + p - MAXLOOP - 2;
927  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
928  minq = MAX2(c0, minq);
929  c0 = j - 3;
930  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
931  maxq = MIN2(c0, maxq);
932  for (q = minq; q < maxq; q++) {
933  if (S[q] != 3)
934  continue;
935 
936  c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
937  ge = MIN2(ge, c0);
938  }
939  }
940  }
941 
942  for (p = i + 2;
943  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
944  p++) {
945  l1 = p - i - 1;
946  if (l1 > MAXLOOP)
947  break;
948 
949  if (S[p] != 3)
950  continue;
951 
952  minq = j - i + p - MAXLOOP - 2;
953  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
954  minq = MAX2(c0, minq);
955  c0 = j - 1;
956  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
957  maxq = MIN2(c0, maxq);
958  for (q = minq; q < maxq; q++) {
959  if (S[q] != 3)
960  continue;
961 
962  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
963  ge = MIN2(ge, c0);
964  }
965  }
966 
967  q = j - 1;
968  if (S[q] == 3)
969  for (p = i + 4;
970  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
971  p++) {
972  l1 = p - i - 1;
973  if (l1 > MAXLOOP)
974  break;
975 
976  if (S[p] != 3)
977  continue;
978 
979  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
980  ge = MIN2(ge, c0);
981  }
982 
983 #if 0
984  /* here comes the additional stuff for the odd dangle models */
985  if (dangles % 1) {
986  en1 = energy + P->dangle5[type][si];
987  en2 = energy + P->dangle5[type][sj];
988  en3 = energy + P->mismatchI[type][si][sj];
989 
990  /* first case with 5' dangle (i.e. j-1) onto enclosing pair */
991  p = i + 1;
992  if (S[p] == 3) {
993  if (p < j - VRNA_GQUAD_MIN_BOX_SIZE) {
994  minq = j - i + p - MAXLOOP - 2;
995  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
996  minq = MAX2(c0, minq);
997  c0 = j - 4;
998  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
999  maxq = MIN2(c0, maxq);
1000  for (q = minq; q < maxq; q++) {
1001  if (S[q] != 3)
1002  continue;
1003 
1004  c0 = en1 + ggg[index[q] + p] + P->internal_loop[j - q - 1];
1005  ge = MIN2(ge, c0);
1006  }
1007  }
1008  }
1009 
1010  for (p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++) {
1011  l1 = p - i - 1;
1012  if (l1 > MAXLOOP)
1013  break;
1014 
1015  if (S[p] != 3)
1016  continue;
1017 
1018  minq = j - i + p - MAXLOOP - 2;
1019  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1020  minq = MAX2(c0, minq);
1021  c0 = j - 2;
1022  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1023  maxq = MIN2(c0, maxq);
1024  for (q = minq; q < maxq; q++) {
1025  if (S[q] != 3)
1026  continue;
1027 
1028  c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
1029  ge = MIN2(ge, c0);
1030  }
1031  }
1032 
1033  q = j - 2;
1034  if (S[q] == 3)
1035  for (p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++) {
1036  l1 = p - i - 1;
1037  if (l1 > MAXLOOP)
1038  break;
1039 
1040  if (S[p] != 3)
1041  continue;
1042 
1043  c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + 1];
1044  ge = MIN2(ge, c0);
1045  }
1046 
1047  /* second case with 3' dangle (i.e. i+1) onto enclosing pair */
1048  }
1049 
1050 #endif
1051  return ge;
1052 }
1053 
1054 
1055 PRIVATE INLINE
1056 int
1057 E_GQuad_IntLoop_comparative(int i,
1058  int j,
1059  unsigned int *tt,
1060  short *S_cons,
1061  short **S5,
1062  short **S3,
1063  unsigned int **a2s,
1064  int *ggg,
1065  int *index,
1066  int n_seq,
1067  vrna_param_t *P)
1068 {
1069  unsigned int type;
1070  int eee, energy, ge, p, q, l1, u1, u2, minq, maxq, c0, s;
1071  vrna_md_t *md;
1072 
1073  md = &(P->model_details);
1074  energy = 0;
1075 
1076  for (s = 0; s < n_seq; s++) {
1077  type = tt[s];
1078  if (md->dangles == 2)
1079  energy += P->mismatchI[type][S3[s][i]][S5[s][j]];
1080 
1081  if (type > 2)
1082  energy += P->TerminalAU;
1083  }
1084 
1085  ge = INF;
1086 
1087  p = i + 1;
1088  if (S_cons[p] == 3) {
1089  if (p < j - VRNA_GQUAD_MIN_BOX_SIZE) {
1090  minq = j - i + p - MAXLOOP - 2;
1091  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1092  minq = MAX2(c0, minq);
1093  c0 = j - 3;
1094  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1095  maxq = MIN2(c0, maxq);
1096  for (q = minq; q < maxq; q++) {
1097  if (S_cons[q] != 3)
1098  continue;
1099 
1100  eee = 0;
1101 
1102  for (s = 0; s < n_seq; s++) {
1103  u1 = a2s[s][j - 1] - a2s[s][q];
1104  eee += P->internal_loop[u1];
1105  }
1106 
1107  c0 = energy +
1108  ggg[index[q] + p] +
1109  eee;
1110  ge = MIN2(ge, c0);
1111  }
1112  }
1113  }
1114 
1115  for (p = i + 2;
1116  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1117  p++) {
1118  l1 = p - i - 1;
1119  if (l1 > MAXLOOP)
1120  break;
1121 
1122  if (S_cons[p] != 3)
1123  continue;
1124 
1125  minq = j - i + p - MAXLOOP - 2;
1126  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1127  minq = MAX2(c0, minq);
1128  c0 = j - 1;
1129  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1130  maxq = MIN2(c0, maxq);
1131  for (q = minq; q < maxq; q++) {
1132  if (S_cons[q] != 3)
1133  continue;
1134 
1135  eee = 0;
1136 
1137  for (s = 0; s < n_seq; s++) {
1138  u1 = a2s[s][p - 1] - a2s[s][i];
1139  u2 = a2s[s][j - 1] - a2s[s][q];
1140  eee += P->internal_loop[u1 + u2];
1141  }
1142 
1143  c0 = energy +
1144  ggg[index[q] + p] +
1145  eee;
1146  ge = MIN2(ge, c0);
1147  }
1148  }
1149 
1150  q = j - 1;
1151  if (S_cons[q] == 3)
1152  for (p = i + 4;
1153  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1154  p++) {
1155  l1 = p - i - 1;
1156  if (l1 > MAXLOOP)
1157  break;
1158 
1159  if (S_cons[p] != 3)
1160  continue;
1161 
1162  eee = 0;
1163 
1164  for (s = 0; s < n_seq; s++) {
1165  u1 = a2s[s][p - 1] - a2s[s][i];
1166  eee += P->internal_loop[u1];
1167  }
1168 
1169  c0 = energy +
1170  ggg[index[q] + p] +
1171  eee;
1172  ge = MIN2(ge, c0);
1173  }
1174 
1175  return ge;
1176 }
1177 
1178 
1179 PRIVATE INLINE
1180 int
1181 E_GQuad_IntLoop_L_comparative(int i,
1182  int j,
1183  unsigned int *tt,
1184  short *S_cons,
1185  short **S5,
1186  short **S3,
1187  unsigned int **a2s,
1188  int **ggg,
1189  int n_seq,
1190  vrna_param_t *P)
1191 {
1192  unsigned int type;
1193  int eee, energy, ge, p, q, l1, u1, u2, minq, maxq, c0, s;
1194  vrna_md_t *md;
1195 
1196  md = &(P->model_details);
1197  energy = 0;
1198 
1199  for (s = 0; s < n_seq; s++) {
1200  type = tt[s];
1201  if (md->dangles == 2)
1202  energy += P->mismatchI[type][S3[s][i]][S5[s][j]];
1203 
1204  if (type > 2)
1205  energy += P->TerminalAU;
1206  }
1207 
1208  ge = INF;
1209 
1210  p = i + 1;
1211  if (S_cons[p] == 3) {
1212  if (p < j - VRNA_GQUAD_MIN_BOX_SIZE) {
1213  minq = j - i + p - MAXLOOP - 2;
1214  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1215  minq = MAX2(c0, minq);
1216  c0 = j - 3;
1217  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1218  maxq = MIN2(c0, maxq);
1219  for (q = minq; q < maxq; q++) {
1220  if (S_cons[q] != 3)
1221  continue;
1222 
1223  eee = 0;
1224 
1225  for (s = 0; s < n_seq; s++) {
1226  u1 = a2s[s][j - 1] - a2s[s][q];
1227  eee += P->internal_loop[u1];
1228  }
1229 
1230  c0 = energy +
1231  ggg[p][q - p] +
1232  eee;
1233  ge = MIN2(ge, c0);
1234  }
1235  }
1236  }
1237 
1238  for (p = i + 2;
1239  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1240  p++) {
1241  l1 = p - i - 1;
1242  if (l1 > MAXLOOP)
1243  break;
1244 
1245  if (S_cons[p] != 3)
1246  continue;
1247 
1248  minq = j - i + p - MAXLOOP - 2;
1249  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1250  minq = MAX2(c0, minq);
1251  c0 = j - 1;
1252  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1253  maxq = MIN2(c0, maxq);
1254  for (q = minq; q < maxq; q++) {
1255  if (S_cons[q] != 3)
1256  continue;
1257 
1258  eee = 0;
1259 
1260  for (s = 0; s < n_seq; s++) {
1261  u1 = a2s[s][p - 1] - a2s[s][i];
1262  u2 = a2s[s][j - 1] - a2s[s][q];
1263  eee += P->internal_loop[u1 + u2];
1264  }
1265 
1266  c0 = energy +
1267  ggg[p][q - p] +
1268  eee;
1269  ge = MIN2(ge, c0);
1270  }
1271  }
1272 
1273  q = j - 1;
1274  if (S_cons[q] == 3)
1275  for (p = i + 4;
1276  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1277  p++) {
1278  l1 = p - i - 1;
1279  if (l1 > MAXLOOP)
1280  break;
1281 
1282  if (S_cons[p] != 3)
1283  continue;
1284 
1285  eee = 0;
1286 
1287  for (s = 0; s < n_seq; s++) {
1288  u1 = a2s[s][p - 1] - a2s[s][i];
1289  eee += P->internal_loop[u1];
1290  }
1291 
1292  c0 = energy +
1293  ggg[p][q - p] +
1294  eee;
1295  ge = MIN2(ge, c0);
1296  }
1297 
1298  return ge;
1299 }
1300 
1301 
1302 PRIVATE INLINE
1303 int *
1304 E_GQuad_IntLoop_exhaustive(int i,
1305  int j,
1306  int **p_p,
1307  int **q_p,
1308  int type,
1309  short *S,
1310  int *ggg,
1311  int threshold,
1312  int *index,
1313  vrna_param_t *P)
1314 {
1315  int energy, *ge, dangles, p, q, l1, minq, maxq, c0;
1316  short si, sj;
1317  int cnt = 0;
1318 
1319  dangles = P->model_details.dangles;
1320  si = S[i + 1];
1321  sj = S[j - 1];
1322  energy = 0;
1323 
1324  if (dangles == 2)
1325  energy += P->mismatchI[type][si][sj];
1326 
1327  if (type > 2)
1328  energy += P->TerminalAU;
1329 
1330  /* guess how many gquads are possible in interval [i+1,j-1] */
1331  *p_p = (int *)vrna_alloc(sizeof(int) * 256);
1332  *q_p = (int *)vrna_alloc(sizeof(int) * 256);
1333  ge = (int *)vrna_alloc(sizeof(int) * 256);
1334 
1335  p = i + 1;
1336  if (S[p] == 3) {
1337  if (p < j - VRNA_GQUAD_MIN_BOX_SIZE) {
1338  minq = j - i + p - MAXLOOP - 2;
1339  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1340  minq = MAX2(c0, minq);
1341  c0 = j - 3;
1342  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1343  maxq = MIN2(c0, maxq);
1344  for (q = minq; q < maxq; q++) {
1345  if (S[q] != 3)
1346  continue;
1347 
1348  c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
1349  if (c0 <= threshold) {
1350  ge[cnt] = energy + P->internal_loop[j - q - 1];
1351  (*p_p)[cnt] = p;
1352  (*q_p)[cnt++] = q;
1353  }
1354  }
1355  }
1356  }
1357 
1358  for (p = i + 2;
1359  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1360  p++) {
1361  l1 = p - i - 1;
1362  if (l1 > MAXLOOP)
1363  break;
1364 
1365  if (S[p] != 3)
1366  continue;
1367 
1368  minq = j - i + p - MAXLOOP - 2;
1369  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1370  minq = MAX2(c0, minq);
1371  c0 = j - 1;
1372  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1373  maxq = MIN2(c0, maxq);
1374  for (q = minq; q < maxq; q++) {
1375  if (S[q] != 3)
1376  continue;
1377 
1378  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
1379  if (c0 <= threshold) {
1380  ge[cnt] = energy + P->internal_loop[l1 + j - q - 1];
1381  (*p_p)[cnt] = p;
1382  (*q_p)[cnt++] = q;
1383  }
1384  }
1385  }
1386 
1387  q = j - 1;
1388  if (S[q] == 3)
1389  for (p = i + 4;
1390  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1391  p++) {
1392  l1 = p - i - 1;
1393  if (l1 > MAXLOOP)
1394  break;
1395 
1396  if (S[p] != 3)
1397  continue;
1398 
1399  c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
1400  if (c0 <= threshold) {
1401  ge[cnt] = energy + P->internal_loop[l1];
1402  (*p_p)[cnt] = p;
1403  (*q_p)[cnt++] = q;
1404  }
1405  }
1406 
1407  (*p_p)[cnt] = -1;
1408 
1409  return ge;
1410 }
1411 
1412 
1413 PRIVATE INLINE
1414 int
1415 E_GQuad_IntLoop_L(int i,
1416  int j,
1417  int type,
1418  short *S,
1419  int **ggg,
1420  int maxdist,
1421  vrna_param_t *P)
1422 {
1423  int energy, ge, dangles, p, q, l1, minq, maxq, c0;
1424  short si, sj;
1425 
1426  dangles = P->model_details.dangles;
1427  si = S[i + 1];
1428  sj = S[j - 1];
1429  energy = 0;
1430 
1431  if (dangles == 2)
1432  energy += P->mismatchI[type][si][sj];
1433 
1434  if (type > 2)
1435  energy += P->TerminalAU;
1436 
1437  ge = INF;
1438 
1439  p = i + 1;
1440  if (S[p] == 3) {
1441  if (p < j - VRNA_GQUAD_MIN_BOX_SIZE) {
1442  minq = j - i + p - MAXLOOP - 2;
1443  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1444  minq = MAX2(c0, minq);
1445  c0 = j - 3;
1446  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1447  maxq = MIN2(c0, maxq);
1448  for (q = minq; q < maxq; q++) {
1449  if (S[q] != 3)
1450  continue;
1451 
1452  c0 = energy + ggg[p][q - p] + P->internal_loop[j - q - 1];
1453  ge = MIN2(ge, c0);
1454  }
1455  }
1456  }
1457 
1458  for (p = i + 2;
1459  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1460  p++) {
1461  l1 = p - i - 1;
1462  if (l1 > MAXLOOP)
1463  break;
1464 
1465  if (S[p] != 3)
1466  continue;
1467 
1468  minq = j - i + p - MAXLOOP - 2;
1469  c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1470  minq = MAX2(c0, minq);
1471  c0 = j - 1;
1472  maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1473  maxq = MIN2(c0, maxq);
1474  for (q = minq; q < maxq; q++) {
1475  if (S[q] != 3)
1476  continue;
1477 
1478  c0 = energy + ggg[p][q - p] + P->internal_loop[l1 + j - q - 1];
1479  ge = MIN2(ge, c0);
1480  }
1481  }
1482 
1483  q = j - 1;
1484  if (S[q] == 3)
1485  for (p = i + 4;
1486  p < j - VRNA_GQUAD_MIN_BOX_SIZE;
1487  p++) {
1488  l1 = p - i - 1;
1489  if (l1 > MAXLOOP)
1490  break;
1491 
1492  if (S[p] != 3)
1493  continue;
1494 
1495  c0 = energy + ggg[p][q - p] + P->internal_loop[l1];
1496  ge = MIN2(ge, c0);
1497  }
1498 
1499  return ge;
1500 }
1501 
1502 
1503 PRIVATE INLINE
1504 FLT_OR_DBL
1505 exp_E_GQuad_IntLoop(int i,
1506  int j,
1507  int type,
1508  short *S,
1509  FLT_OR_DBL *G,
1510  FLT_OR_DBL *scale,
1511  int *index,
1512  vrna_exp_param_t *pf)
1513 {
1514  int k, l, minl, maxl, u, r;
1515  FLT_OR_DBL q, qe;
1516  double *expintern;
1517  short si, sj;
1518 
1519  q = 0;
1520  si = S[i + 1];
1521  sj = S[j - 1];
1522  qe = (FLT_OR_DBL)pf->expmismatchI[type][si][sj];
1523  expintern = &(pf->expinternal[0]);
1524 
1525  if (type > 2)
1526  qe *= (FLT_OR_DBL)pf->expTermAU;
1527 
1528  k = i + 1;
1529  if (S[k] == 3) {
1530  if (k < j - VRNA_GQUAD_MIN_BOX_SIZE) {
1531  minl = j - MAXLOOP - 1;
1532  u = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1533  minl = MAX2(u, minl);
1534  u = j - 3;
1535  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1536  maxl = MIN2(u, maxl);
1537  for (l = minl; l < maxl; l++) {
1538  if (S[l] != 3)
1539  continue;
1540 
1541  if (G[index[k] - l] == 0.)
1542  continue;
1543 
1544  q += qe
1545  * G[index[k] - l]
1546  * (FLT_OR_DBL)expintern[j - l - 1]
1547  * scale[j - l + 1];
1548  }
1549  }
1550  }
1551 
1552  for (k = i + 2;
1553  k <= j - VRNA_GQUAD_MIN_BOX_SIZE;
1554  k++) {
1555  u = k - i - 1;
1556  if (u > MAXLOOP)
1557  break;
1558 
1559  if (S[k] != 3)
1560  continue;
1561 
1562  minl = j - i + k - MAXLOOP - 2;
1563  r = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1564  minl = MAX2(r, minl);
1565  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1566  r = j - 1;
1567  maxl = MIN2(r, maxl);
1568  for (l = minl; l < maxl; l++) {
1569  if (S[l] != 3)
1570  continue;
1571 
1572  if (G[index[k] - l] == 0.)
1573  continue;
1574 
1575  q += qe
1576  * G[index[k] - l]
1577  * (FLT_OR_DBL)expintern[u + j - l - 1]
1578  * scale[u + j - l + 1];
1579  }
1580  }
1581 
1582  l = j - 1;
1583  if (S[l] == 3)
1584  for (k = i + 4; k <= j - VRNA_GQUAD_MIN_BOX_SIZE; k++) {
1585  u = k - i - 1;
1586  if (u > MAXLOOP)
1587  break;
1588 
1589  if (S[k] != 3)
1590  continue;
1591 
1592  if (G[index[k] - l] == 0.)
1593  continue;
1594 
1595  q += qe
1596  * G[index[k] - l]
1597  * (FLT_OR_DBL)expintern[u]
1598  * scale[u + 2];
1599  }
1600 
1601  return q;
1602 }
1603 
1604 
1605 PRIVATE INLINE
1606 FLT_OR_DBL
1607 exp_E_GQuad_IntLoop_comparative(int i,
1608  int j,
1609  unsigned int *tt,
1610  short *S_cons,
1611  short **S5,
1612  short **S3,
1613  unsigned int **a2s,
1614  FLT_OR_DBL *G,
1615  FLT_OR_DBL *scale,
1616  int *index,
1617  int n_seq,
1618  vrna_exp_param_t *pf)
1619 {
1620  unsigned int type;
1621  int k, l, minl, maxl, u, u1, u2, r, s;
1622  FLT_OR_DBL q, qe, qqq;
1623  double *expintern;
1624  vrna_md_t *md;
1625 
1626  q = 0;
1627  qe = 1.;
1628  md = &(pf->model_details);
1629  expintern = &(pf->expinternal[0]);
1630 
1631  for (s = 0; s < n_seq; s++) {
1632  type = tt[s];
1633  if (md->dangles == 2)
1634  qe *= (FLT_OR_DBL)pf->expmismatchI[type][S3[s][i]][S5[s][j]];
1635 
1636  if (type > 2)
1637  qe *= (FLT_OR_DBL)pf->expTermAU;
1638  }
1639 
1640  k = i + 1;
1641  if (S_cons[k] == 3) {
1642  if (k < j - VRNA_GQUAD_MIN_BOX_SIZE) {
1643  minl = j - MAXLOOP - 1;
1644  u = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1645  minl = MAX2(u, minl);
1646  u = j - 3;
1647  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1648  maxl = MIN2(u, maxl);
1649  for (l = minl; l < maxl; l++) {
1650  if (S_cons[l] != 3)
1651  continue;
1652 
1653  if (G[index[k] - l] == 0.)
1654  continue;
1655 
1656  qqq = 1.;
1657 
1658  for (s = 0; s < n_seq; s++) {
1659  u1 = a2s[s][j - 1] - a2s[s][l];
1660  qqq *= expintern[u1];
1661  }
1662 
1663  q += qe *
1664  G[index[k] - l] *
1665  qqq *
1666  scale[j - l + 1];
1667  }
1668  }
1669  }
1670 
1671  for (k = i + 2;
1672  k <= j - VRNA_GQUAD_MIN_BOX_SIZE;
1673  k++) {
1674  u = k - i - 1;
1675  if (u > MAXLOOP)
1676  break;
1677 
1678  if (S_cons[k] != 3)
1679  continue;
1680 
1681  minl = j - i + k - MAXLOOP - 2;
1682  r = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
1683  minl = MAX2(r, minl);
1684  maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
1685  r = j - 1;
1686  maxl = MIN2(r, maxl);
1687  for (l = minl; l < maxl; l++) {
1688  if (S_cons[l] != 3)
1689  continue;
1690 
1691  if (G[index[k] - l] == 0.)
1692  continue;
1693 
1694  qqq = 1.;
1695 
1696  for (s = 0; s < n_seq; s++) {
1697  u1 = a2s[s][k - 1] - a2s[s][i];
1698  u2 = a2s[s][j - 1] - a2s[s][l];
1699  qqq *= expintern[u1 + u2];
1700  }
1701 
1702  q += qe *
1703  G[index[k] - l] *
1704  qqq *
1705  scale[u + j - l + 1];
1706  }
1707  }
1708 
1709  l = j - 1;
1710  if (S_cons[l] == 3)
1711  for (k = i + 4; k <= j - VRNA_GQUAD_MIN_BOX_SIZE; k++) {
1712  u = k - i - 1;
1713  if (u > MAXLOOP)
1714  break;
1715 
1716  if (S_cons[k] != 3)
1717  continue;
1718 
1719  if (G[index[k] - l] == 0.)
1720  continue;
1721 
1722  qqq = 1.;
1723 
1724  for (s = 0; s < n_seq; s++) {
1725  u1 = a2s[s][k - 1] - a2s[s][i];
1726  qqq *= expintern[u1];
1727  }
1728 
1729  q += qe *
1730  G[index[k] - l] *
1731  qqq *
1732  scale[u + 2];
1733  }
1734 
1735  return q;
1736 }
1737 
1738 
1744 #endif
vrna_alloc
void * vrna_alloc(unsigned size)
Allocate space safely.
vrna_fc_s::S
short ** S
Numerical encoding of the sequences in the alignment.
Definition: fold_compound.h:271
vrna_fc_s::S_cons
short * S_cons
Numerical encoding of the consensus sequence.
Definition: fold_compound.h:268
vrna_fc_s::ptype
char * ptype
Pair type array.
Definition: fold_compound.h:226
vrna_fc_s::sequence_encoding
short * sequence_encoding
Numerical encoding of the input sequence.
Definition: fold_compound.h:221
vrna_param_s::model_details
vrna_md_t model_details
Model details to be used in the recursions.
Definition: basic.h:96
vrna_exp_param_s::model_details
vrna_md_t model_details
Model details to be used in the recursions.
Definition: basic.h:154
FLT_OR_DBL
double FLT_OR_DBL
Typename for floating point number in partition function computations.
Definition: basic.h:43
vrna_fc_s
The most basic data structure required by many functions throughout the RNAlib.
Definition: fold_compound.h:132
vrna_param_s
The datastructure that contains temperature scaled energy parameters.
Definition: basic.h:57
MIN2
#define MIN2(A, B)
Get the minimum of two comparable values.
Definition: basic.h:111
MAXLOOP
#define MAXLOOP
Definition: constants.h:29
parse_gquad
int parse_gquad(const char *struc, int *L, int l[3])
INF
#define INF
Definition: constants.h:17
vrna_fc_s::jindx
int * jindx
DP matrix accessor.
Definition: fold_compound.h:167
basic.h
Various data structures and pre-processor macros.
backtrack_GQuad_IntLoop
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:423
vrna_md_s
The data structure that contains the complete model details used throughout the calculations.
Definition: model.h:180
vrna_exp_param_s
The data structure that contains temperature scaled Boltzmann weights of the energy parameters...
Definition: basic.h:103
VRNA_FC_TYPE_COMPARATIVE
Definition: fold_compound.h:116
MAX2
#define MAX2(A, B)
Get the maximum of two comparable values.
Definition: basic.h:116
vrna_fc_s::params
vrna_param_t * params
The precomputed free energy contributions for each type of loop.
Definition: fold_compound.h:163
vrna_fc_s::n_seq
unsigned int n_seq
The number of sequences in the alignment.
Definition: fold_compound.h:262
vrna_fc_s::matrices
vrna_mx_mfe_t * matrices
The MFE DP matrices.
Definition: fold_compound.h:160
VRNA_FC_TYPE_SINGLE
Definition: fold_compound.h:115
vrna_md_s::dangles
int dangles
Specifies the dangle model used in any energy evaluation (0,1,2 or 3)
Definition: model.h:184
vrna_bp_stack_s
Base pair stack element.
Definition: basic.h:143
pr
FLT_OR_DBL * pr
A pointer to the base pair probability matrix.
basic.h
Functions to deal with sets of energy parameters.
fold_compound.h
The Basic Fold Compound API.
dangles
int dangles
Switch the energy model for dangling end contributions (0, 1, 2, 3)
get_gquad_matrix
int * get_gquad_matrix(short *S, vrna_param_t *P)
Get a triangular matrix prefilled with minimum free energy contributions of G-quadruplexes.
vrna_elem_prob_s
Data structure representing a single entry of an element probability list (e.g. list of pair probabil...
Definition: structures.h:453
vrna_mx_mfe_s::ggg
int * ggg
Energies of g-quadruplexes.
Definition: dp_matrices.h:73
vrna_fc_s::type
const vrna_fc_type_e type
The type of the vrna_fold_compound_t.
Definition: fold_compound.h:137
backtrack_GQuad_IntLoop_L
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:666