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
16 /**
17 * @file gquad.h
18 * @ingroup paired_modules
19 * @brief G-quadruplexes
20 */
21
22 /**
23 * @addtogroup gquads
24 * @{
25 *
26 * @brief Various functions related to G-quadruplex computations
27 */
28
29
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
50 /**
51 * @brief Get a triangular matrix prefilled with minimum free energy
52 * contributions of G-quadruplexes.
53 *
54 * At each position ij in the matrix, the minimum free energy of any
55 * G-quadruplex delimited by i and j is stored. If no G-quadruplex formation
56 * is possible, the matrix element is set to INF.
57 * Access the elements in the matrix via matrix[indx[j]+i]. To get
58 * the integer array indx see get_jindx().
59 *
60 * @see get_jindx(), encode_sequence()
61 *
62 * @param S The encoded sequence
63 * @param P A pointer to the data structure containing the precomputed energy contributions
64 * @return A pointer to the G-quadruplex contribution matrix
65 */
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
189 /**
190 * given a dot-bracket structure (possibly) containing gquads encoded
191 * by '+' signs, find first gquad, return end position or 0 if none found
192 * Upon return L and l[] contain the number of stacked layers, as well as
193 * the lengths of the linker regions.
194 * To parse a string with many gquads, call parse_gquad repeatedly e.g.
195 * end1 = parse_gquad(struc, &L, l); ... ;
196 * end2 = parse_gquad(struc+end1, &L, l); end2+=end1; ... ;
197 * end3 = parse_gquad(struc+end2, &L, l); end3+=end2; ... ;
198 */
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
vrna_BT_gquad_mfe(vrna_fold_compound_t * vc,int i,int j,vrna_bp_stack_t * bp_stack,int * stack_count)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
278 case VRNA_FC_TYPE_COMPARATIVE:
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
vrna_BT_gquad_int(vrna_fold_compound_t * vc,int i,int j,int en,vrna_bp_stack_t * bp_stack,int * stack_count)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
407 /**
408 * backtrack an interior loop like enclosed g-quadruplex
409 * with closing pair (i,j)
410 *
411 * @param c The total contribution the loop should resemble
412 * @param i position i of enclosing pair
413 * @param j position j of enclosing pair
414 * @param type base pair type of enclosing pair (must be reverse type)
415 * @param S integer encoded sequence
416 * @param ggg triangular matrix containing g-quadruplex contributions
417 * @param index the index for accessing the triangular matrix
418 * @param p here the 5' position of the gquad is stored
419 * @param q here the 3' position of the gquad is stored
420 * @param P the datastructure containing the precalculated contibutions
421 *
422 * @return 1 on success, 0 if no gquad found
423 */
424 INLINE 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)425 backtrack_GQuad_IntLoop(int c,
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
439 dangles = P->model_details.dangles;
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
backtrack_GQuad_IntLoop_comparative(int c,int i,int j,unsigned int * type,short * S_cons,short ** S5,short ** S3,unsigned int ** a2s,int * ggg,int * index,int * p,int * q,int n_seq,vrna_param_t * P)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
541 dangles = P->model_details.dangles;
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
651 /**
652 * backtrack an interior loop like enclosed g-quadruplex
653 * with closing pair (i,j) with underlying Lfold matrix
654 *
655 * @param c The total contribution the loop should resemble
656 * @param i position i of enclosing pair
657 * @param j position j of enclosing pair
658 * @param type base pair type of enclosing pair (must be reverse type)
659 * @param S integer encoded sequence
660 * @param ggg triangular matrix containing g-quadruplex contributions
661 * @param p here the 5' position of the gquad is stored
662 * @param q here the 3' position of the gquad is stored
663 * @param P the datastructure containing the precalculated contibutions
664 *
665 * @return 1 on success, 0 if no gquad found
666 */
667 INLINE 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)668 backtrack_GQuad_IntLoop_L(int c,
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
682 dangles = P->model_details.dangles;
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
backtrack_GQuad_IntLoop_L_comparative(int c,int i,int j,unsigned int * type,short * S_cons,short ** S5,short ** S3,unsigned int ** a2s,int ** ggg,int * p,int * q,int n_seq,vrna_param_t * P)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
E_GQuad_IntLoop(int i,int j,int type,short * S,int * ggg,int * index,vrna_param_t * P)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
912 dangles = P->model_details.dangles;
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
E_GQuad_IntLoop_comparative(int i,int j,unsigned int * tt,short * S_cons,short ** S5,short ** S3,unsigned int ** a2s,int * ggg,int * index,int n_seq,vrna_param_t * P)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
E_GQuad_IntLoop_L_comparative(int i,int j,unsigned int * tt,short * S_cons,short ** S5,short ** S3,unsigned int ** a2s,int ** ggg,int n_seq,vrna_param_t * P)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 *
E_GQuad_IntLoop_exhaustive(int i,int j,int ** p_p,int ** q_p,int type,short * S,int * ggg,int threshold,int * index,vrna_param_t * P)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
1321 dangles = P->model_details.dangles;
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
E_GQuad_IntLoop_L(int i,int j,int type,short * S,int ** ggg,int maxdist,vrna_param_t * P)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
1428 dangles = P->model_details.dangles;
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
exp_E_GQuad_IntLoop(int i,int j,int type,short * S,FLT_OR_DBL * G,FLT_OR_DBL * scale,int * index,vrna_exp_param_t * pf)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
exp_E_GQuad_IntLoop_comparative(int i,int j,unsigned int * tt,short * S_cons,short ** S5,short ** S3,unsigned int ** a2s,FLT_OR_DBL * G,FLT_OR_DBL * scale,int * index,int n_seq,vrna_exp_param_t * pf)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
1741 /**
1742 * @}
1743 */
1744
1745
1746 #endif
1747