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