1 /*
2 * compute potentially pseudoknotted structures of an RNA sequence
3 * Ivo Hofacker
4 * Vienna RNA package
5 */
6
7 /*
8 * library containing the function used in PKplex
9 * it generates pseudoknotted structures by letting the sequence form a duplex structure with itself
10 */
11
12 #ifdef HAVE_CONFIG_H
13 #include "config.h"
14 #endif
15
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <math.h>
19 #include <ctype.h>
20 #include <string.h>
21 #include <time.h>
22
23 #include "ViennaRNA/params/default.h"
24 #include "ViennaRNA/fold_vars.h"
25 #include "ViennaRNA/utils/basic.h"
26 #include "ViennaRNA/params/basic.h"
27 #include "ViennaRNA/loops/all.h"
28 #include "ViennaRNA/alphabet.h"
29 #include "ViennaRNA/LPfold.h"
30 #include "ViennaRNA/mfe.h"
31 #include "ViennaRNA/datastructures/heap.h"
32
33 #include "ViennaRNA/loops/external_hc.inc"
34
35 #include "ViennaRNA/pk_plex.h"
36
37 #undef MAXLOOP
38 #define MAXLOOP 10
39
40 #define DEFAULT_PENALTY 810
41 #define DEFAULT_DELTA 0
42 #define DEFAULT_INTERACTION_LENGTH 12
43 #define DEFAULT_CUTOFF 1e-3
44
45 struct vrna_pk_plex_option_s {
46 unsigned int delta;
47 unsigned int max_interaction_length;
48 int min_penalty;
49 vrna_callback_pk_plex_score *scoring_function;
50 void *scoring_data;
51 };
52
53 typedef struct {
54 int penalty;
55 } default_data;
56
57 /*
58 #################################
59 # GLOBAL VARIABLES #
60 #################################
61 */
62
63 /*
64 #################################
65 # PRIVATE VARIABLES #
66 #################################
67 */
68
69 /*
70 #################################
71 # PRIVATE FUNCTION DECLARATIONS #
72 #################################
73 */
74 PRIVATE int
75 default_pk_plex_penalty(const short *pt,
76 int i,
77 int j,
78 int k,
79 int l,
80 void *data);
81
82
83 PRIVATE vrna_heap_t
84 duplexfold_XS(vrna_fold_compound_t *fc,
85 const int **access_s1,
86 const int max_interaction_length,
87 vrna_callback_pk_plex_score *scoring_function,
88 void *scoring_data);
89
90
91 PRIVATE char *
92 backtrack_XS(vrna_fold_compound_t *fc,
93 int kk,
94 int ll,
95 const int ii,
96 const int jj,
97 const int max_interaction_length,
98 int ***c3);
99
100
101 PRIVATE int
102 prepare_PKplex(vrna_fold_compound_t *fc);
103
104
105 PRIVATE int
106 PKplex_heap_cmp(const void *a,
107 const void *b,
108 void *data);
109
110
111 /*
112 #################################
113 # BEGIN OF FUNCTION DEFINITIONS #
114 #################################
115 */
116 PUBLIC vrna_pk_plex_opt_t
vrna_pk_plex_opt_defaults(void)117 vrna_pk_plex_opt_defaults(void)
118 {
119 return vrna_pk_plex_opt(DEFAULT_DELTA,
120 DEFAULT_INTERACTION_LENGTH,
121 DEFAULT_PENALTY);
122 }
123
124
125 PUBLIC vrna_pk_plex_opt_t
vrna_pk_plex_opt(unsigned int delta,unsigned int max_interaction_length,int pk_penalty)126 vrna_pk_plex_opt(unsigned int delta,
127 unsigned int max_interaction_length,
128 int pk_penalty)
129 {
130 vrna_pk_plex_opt_t opt;
131
132 opt = (vrna_pk_plex_opt_t)vrna_alloc(sizeof(struct vrna_pk_plex_option_s));
133
134 opt->delta = delta;
135 opt->max_interaction_length = max_interaction_length;
136 opt->min_penalty = pk_penalty;
137 opt->scoring_function = NULL;
138 opt->scoring_data = NULL;
139
140 return opt;
141 }
142
143
144 PUBLIC vrna_pk_plex_opt_t
vrna_pk_plex_opt_fun(unsigned int delta,unsigned int max_interaction_length,vrna_callback_pk_plex_score * scoring_function,void * scoring_data)145 vrna_pk_plex_opt_fun(unsigned int delta,
146 unsigned int max_interaction_length,
147 vrna_callback_pk_plex_score *scoring_function,
148 void *scoring_data)
149 {
150 vrna_pk_plex_opt_t opt = NULL;
151
152 if (scoring_function) {
153 opt = (vrna_pk_plex_opt_t)vrna_alloc(sizeof(struct vrna_pk_plex_option_s));
154
155 opt->delta = delta;
156 opt->max_interaction_length = max_interaction_length;
157 opt->scoring_function = scoring_function;
158 opt->scoring_data = scoring_data;
159 }
160
161 return opt;
162 }
163
164
165 PUBLIC vrna_pk_plex_t *
vrna_pk_plex(vrna_fold_compound_t * fc,const int ** accessibility,vrna_pk_plex_opt_t user_options)166 vrna_pk_plex(vrna_fold_compound_t *fc,
167 const int **accessibility,
168 vrna_pk_plex_opt_t user_options)
169 {
170 size_t cnt;
171 char *mfe_struct, *constraint;
172 int i, **opening_energies;
173 double mfe, mfe_pk, constrainedEnergy, pk_penalty, subopts, best_e;
174 default_data scoring_dat;
175 vrna_pk_plex_t *hits, *hit_ptr, *ptr, *mfe_entry;
176 vrna_pk_plex_opt_t options;
177 vrna_heap_t interactions, results;
178
179 hits = NULL;
180 results = NULL;
181 opening_energies = NULL;
182
183 if (fc) {
184 mfe_struct = (char *)vrna_alloc(sizeof(char) * (fc->length + 1));
185
186 mfe = mfe_pk = (double)vrna_mfe(fc, mfe_struct);
187
188 options = (user_options) ? user_options : vrna_pk_plex_opt_defaults();
189
190 /* apply simplified scoring function with constant penalty */
191 if (!options->scoring_function) {
192 scoring_dat.penalty = options->min_penalty;
193 options->scoring_function = &default_pk_plex_penalty;
194 options->scoring_data = (void *)&scoring_dat;
195 }
196
197 /* sanity check for maximum interaction length */
198 options->max_interaction_length = MIN2(DEFAULT_INTERACTION_LENGTH, fc->length - 3);
199
200 /* compute opening energies if not passed as argument */
201 if (!accessibility)
202 opening_energies = vrna_pk_plex_accessibility(fc->sequence,
203 options->max_interaction_length,
204 DEFAULT_CUTOFF);
205
206 /* obtain list of potential PK interactions */
207 interactions = duplexfold_XS(fc,
208 (accessibility) ? accessibility : (const int **)opening_energies,
209 options->max_interaction_length,
210 options->scoring_function,
211 options->scoring_data);
212
213 pk_penalty =
214 ((double)options->scoring_function(NULL, 0, 0, 0, 0, options->scoring_data) / 100.);
215 subopts = ((double)options->delta / 100.);
216
217 if (vrna_heap_size(interactions) > 0) {
218 results = vrna_heap_init(vrna_heap_size(interactions) + 2,
219 PKplex_heap_cmp,
220 NULL,
221 NULL,
222 NULL);
223
224 /* now we re-evaluate the energies and thereby prune the list of pre-results
225 * such that the re-avaluation process is not done too often.
226 */
227
228 constraint = (char *)vrna_alloc(sizeof(char) * (fc->length + 1));
229
230 while ((hit_ptr = vrna_heap_pop(interactions))) {
231 /*
232 * simple check whether we believe that this structure might achieve
233 * a net free energy within our boundaries
234 * For that, we add up the interaction energy, the pk-free MFE
235 * the PK penalty, and the minimal energy to unfold at least
236 * one of the PK interaction sites. This for now seems a good
237 * choice for a lower bound of the actual energy
238 */
239 best_e = hit_ptr->dGint +
240 mfe +
241 pk_penalty +
242 MIN2(hit_ptr->dG1, hit_ptr->dG2);
243
244 if (best_e <= mfe_pk + subopts) {
245 /* now for the exact evaluation of the structures energy incl. PKs */
246
247 /* prepare the structure constraint for constrained folding */
248 vrna_hc_init(fc);
249
250 for (i = hit_ptr->start_5; i <= hit_ptr->end_5; i++)
251 vrna_hc_add_up(fc, i, VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS);
252 for (i = hit_ptr->start_3; i <= hit_ptr->end_3; i++)
253 vrna_hc_add_up(fc, i, VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS);
254
255 /* energy evaluation */
256 constrainedEnergy = vrna_mfe(fc, constraint);
257
258 if (options->scoring_function != &default_pk_plex_penalty) {
259 short *pt = vrna_ptable(constraint);
260 hit_ptr->dGpk = ((double)options->scoring_function((const short *)pt,
261 hit_ptr->start_5,
262 hit_ptr->end_5,
263 hit_ptr->start_3,
264 hit_ptr->end_3,
265 options->scoring_data) / 100.);
266 free(pt);
267 } else {
268 hit_ptr->dGpk = pk_penalty;
269 }
270
271 /*
272 * compute net free energy, i.e.:
273 * interaction energy +
274 * constrained MFE +
275 * pseudo-knot penalty
276 */
277 hit_ptr->energy = hit_ptr->dGint +
278 constrainedEnergy +
279 hit_ptr->dGpk;
280
281 /* check if this structure is worth keeping */
282 if (hit_ptr->energy <= mfe_pk + subopts) {
283 /* add pseudo-knot brackets to the structure */
284 for (i = hit_ptr->start_5 - 1; i < hit_ptr->end_5; i++)
285 if (hit_ptr->structure[i - hit_ptr->start_5 + 1] == '(')
286 constraint[i] = '[';
287
288 for (i = hit_ptr->start_3 - 1; i < hit_ptr->end_3; i++)
289 if (hit_ptr->structure[i - hit_ptr->start_3 + 1 + 1 + 1 +
290 hit_ptr->end_5 - hit_ptr->start_5] == ')')
291 constraint[i] = ']';
292
293 if (hit_ptr->energy < mfe_pk)
294 mfe_pk = hit_ptr->energy;
295
296 free(hit_ptr->structure);
297 hit_ptr->structure = constraint;
298
299 vrna_heap_insert(results, hit_ptr);
300
301 constraint = (char *)vrna_alloc(sizeof(char) * (fc->length + 1));
302
303 continue;
304 }
305 }
306
307 free(hit_ptr->structure);
308 free(hit_ptr);
309 }
310
311 free(constraint);
312 }
313
314 /* Add the PK-free MFE as potential result */
315 mfe_entry = (vrna_pk_plex_t *)vrna_alloc(sizeof(vrna_pk_plex_t));
316
317 mfe_entry->structure = mfe_struct;
318 mfe_entry->energy = mfe;
319 mfe_entry->start_5 = 0;
320
321 if (results == NULL) {
322 results = vrna_heap_init(1,
323 PKplex_heap_cmp,
324 NULL,
325 NULL,
326 NULL);
327 }
328
329 vrna_heap_insert(results, (void *)mfe_entry);
330
331 /*
332 * now go through the active hits again and filter those out that are above
333 * the subopt threshold. This is necessary due to the fact that while we've
334 * been processing the list once we always compared against the best known
335 * pk-MFE. But this value might have changed during processing.
336 */
337 cnt = 0;
338 hits = (vrna_pk_plex_t *)vrna_alloc(sizeof(vrna_pk_plex_t) * (vrna_heap_size(results) + 1));
339 /* collect all final results */
340 while ((ptr = vrna_heap_pop(results))) {
341 if (ptr->energy > mfe_pk + subopts)
342 break;
343
344 hits[cnt++] = *ptr;
345 }
346
347 /* end of list marker */
348 hits[cnt].structure = NULL;
349
350 /* remove intermediate hits that didn't surpass the threshold */
351 while ((ptr = vrna_heap_pop(results))) {
352 free(ptr->structure);
353 free(ptr);
354 }
355
356 /* cleanup the heap storages */
357 vrna_heap_free(interactions);
358 vrna_heap_free(results);
359
360 if (opening_energies) {
361 i = opening_energies[0][0];
362 while (--i > -1)
363 free(opening_energies[i]);
364 free(opening_energies);
365 }
366
367 if (options != user_options)
368 free(options);
369 }
370
371 return hits;
372 }
373
374
375 PUBLIC int **
vrna_pk_plex_accessibility(const char * sequence,unsigned int unpaired,double cutoff)376 vrna_pk_plex_accessibility(const char *sequence,
377 unsigned int unpaired,
378 double cutoff)
379 {
380 unsigned int n, i, j;
381 int **a = NULL;
382 double **pup, kT;
383 plist *dpp = NULL;
384 vrna_fold_compound_t *fc;
385 vrna_param_t *P;
386 vrna_md_t *md;
387
388 if (sequence) {
389 fc = vrna_fold_compound(sequence, NULL, VRNA_OPTION_DEFAULT | VRNA_OPTION_WINDOW);
390
391 n = fc->length;
392 P = fc->params;
393 md = &(P->model_details);
394
395 pup = (double **)vrna_alloc((n + 1) * sizeof(double *));
396 pup[0] = (double *)vrna_alloc(sizeof(double)); /*I only need entry 0*/
397 pup[0][0] = (double)unpaired;
398
399 (void)pfl_fold(fc->sequence, n, n, cutoff, pup, &dpp, NULL, NULL);
400
401 kT = (md->temperature + K0) * GASCONST / 1000.0;
402
403 /* prepare the accesibility array */
404 a = (int **)vrna_alloc(sizeof(int *) * (unpaired + 2));
405
406 for (i = 0; i < unpaired + 2; i++)
407 a[i] = (int *)vrna_alloc(sizeof(int) * (n + 1));
408
409 for (i = 0; i <= n; i++)
410 for (j = 0; j < unpaired + 2; j++)
411 a[j][i] = INF;
412
413 for (i = 1; i <= n; i++) {
414 for (j = 1; j < unpaired + 1; j++)
415 if (pup[i][j] > 0)
416 a[j][i] = rint(100 * (-log(pup[i][j])) * kT);
417 }
418
419 a[0][0] = unpaired + 2;
420
421 vrna_fold_compound_free(fc);
422
423 for (i = 0; i <= n; i++)
424 free(pup[i]);
425 free(pup);
426 }
427
428 return a;
429 }
430
431
432 /*
433 #####################################
434 # BEGIN OF STATIC HELPER FUNCTIONS #
435 #####################################
436 */
437 PRIVATE int
default_pk_plex_penalty(const short * pt,int i,int j,int k,int l,void * data)438 default_pk_plex_penalty(const short *pt,
439 int i,
440 int j,
441 int k,
442 int l,
443 void *data)
444 {
445 return ((default_data *)data)->penalty;
446 }
447
448
449 PRIVATE int ***
get_array(unsigned int n,unsigned int interaction_length)450 get_array(unsigned int n,
451 unsigned int interaction_length)
452 {
453 unsigned int i, j;
454 int ***c3;
455
456 c3 = (int ***)vrna_alloc(sizeof(int **) * (n));
457
458 for (i = 0; i < n; i++) {
459 c3[i] = (int **)vrna_alloc(sizeof(int *) * interaction_length);
460
461 for (j = 0; j < interaction_length; j++)
462 c3[i][j] = (int *)vrna_alloc(sizeof(int) * interaction_length);
463 }
464
465 return c3;
466 }
467
468
469 PRIVATE void
reset_array(int *** c3,unsigned int n,unsigned int interaction_length)470 reset_array(int ***c3,
471 unsigned int n,
472 unsigned int interaction_length)
473 {
474 unsigned int i, j, k;
475
476 for (i = 0; i < n; i++) {
477 for (j = 0; j < interaction_length; j++)
478 for (k = 0; k < interaction_length; k++)
479 c3[i][j][k] = INF;
480 }
481 }
482
483
484 PRIVATE void
free_array(int *** c3,unsigned int n,unsigned int interaction_length)485 free_array(int ***c3,
486 unsigned int n,
487 unsigned int interaction_length)
488 {
489 unsigned int i, j;
490
491 for (i = 0; i < n; i++) {
492 for (j = 0; j < interaction_length; j++)
493 free(c3[i][j]);
494 free(c3[i]);
495 }
496
497 free(c3);
498 }
499
500
501 PRIVATE int
prepare_PKplex(vrna_fold_compound_t * fc)502 prepare_PKplex(vrna_fold_compound_t *fc)
503 {
504 /* prepare Boltzmann factors if required */
505 vrna_params_prepare(fc, VRNA_OPTION_MFE);
506
507 /* prepare ptype array(s) */
508 vrna_ptypes_prepare(fc, VRNA_OPTION_MFE);
509
510 /* prepare hard constraints */
511 vrna_hc_prepare(fc, VRNA_OPTION_MFE);
512
513 /* prepare soft constraints data structure, if required */
514 vrna_sc_prepare(fc, VRNA_OPTION_MFE);
515
516 return 1;
517 }
518
519
520 PRIVATE vrna_heap_t
duplexfold_XS(vrna_fold_compound_t * fc,const int ** access_s1,const int max_interaction_length,vrna_callback_pk_plex_score * scoring_function,void * scoring_data)521 duplexfold_XS(vrna_fold_compound_t *fc,
522 const int **access_s1,
523 const int max_interaction_length,
524 vrna_callback_pk_plex_score *scoring_function,
525 void *scoring_data)
526 {
527 char *struc;
528 short *S, *S1, si, sk, sl, sp, sq;
529 unsigned int n, type, type2, type3;
530 int ***c3, i, j, k, l, p, q, Emin, l_min, k_min, j_min, E,
531 tempK, i_pos_begin, j_pos_end, dGx, dGy, inter,
532 turn, penalty;
533
534 vrna_param_t *P;
535 vrna_md_t *md;
536 vrna_hc_t *hc;
537 vrna_heap_t storage;
538 vrna_pk_plex_t *entry;
539
540 vrna_callback_hc_evaluate *evaluate_ext;
541 struct hc_ext_def_dat hc_dat_local;
542
543 struc = NULL;
544 n = fc->length;
545 S = fc->sequence_encoding2;
546 S1 = fc->sequence_encoding;
547 P = fc->params;
548 md = &(P->model_details);
549 turn = md->min_loop_size;
550 hc = fc->hc;
551 penalty = scoring_function(NULL, 0, 0, 0, 0, scoring_data);
552
553 storage = vrna_heap_init(128,
554 PKplex_heap_cmp,
555 NULL,
556 NULL,
557 NULL);
558
559 evaluate_ext = prepare_hc_ext_def(fc, &hc_dat_local);
560
561 c3 = get_array(n, max_interaction_length);
562
563 if (n > turn + 1) {
564 for (i = n - turn - 1; i > 0; i--) {
565 Emin = INF;
566 j_min = 0;
567 l_min = 0;
568 k_min = 0;
569
570 reset_array(c3, n, max_interaction_length);
571
572 si = S1[i + 1];
573
574 /* matrix starting values for (i,j)-basepairs */
575 for (j = i + turn + 1; j <= n; j++) {
576 if (evaluate_ext(i, j, i, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
577 type = md->pair[S[j]][S[i]];
578 c3[j - 1][max_interaction_length - 1][0] = vrna_E_ext_stem(type,
579 S1[j - 1],
580 si,
581 P);
582 }
583 }
584
585 i_pos_begin = MAX2(0, i - max_interaction_length);
586
587 /* fill matrix */
588 for (k = i - 1; k > i_pos_begin; k--) {
589 tempK = max_interaction_length - i + k - 1;
590 sk = S1[k + 1];
591 for (l = i + turn + 1; l <= n; l++) {
592 if (hc->mx[n * k + l] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) {
593 type2 = md->pair[S[k]][S[l]];
594 sl = S1[l - 1];
595
596 for (p = k + 1; (p <= i) && (p <= k + MAXLOOP + 1); p++) {
597 sp = S1[p - 1];
598
599 for (q = l - 1; (q >= i + turn + 1) && (q >= l - MAXLOOP - 1); q--) {
600 if (p - k + l - q - 2 > MAXLOOP)
601 break;
602
603 if (hc->mx[n * p + q] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC) {
604 type3 = md->pair[S[q]][S[p]];
605 sq = S1[q + 1];
606
607 E = E_IntLoop(p - k - 1,
608 l - q - 1,
609 type2,
610 type3,
611 sk,
612 sl,
613 sp,
614 sq,
615 P);
616 for (j = MAX2(i + turn + 1, l - max_interaction_length + 1); j <= q; j++) {
617 if (hc->mx[n * i + j] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
618 type = md->pair[S[i]][S[j]];
619 c3[j - 1][tempK][l - j] =
620 MIN2(c3[j - 1][tempK][l - j],
621 c3[j - 1][max_interaction_length - i + p - 1][q - j] + E);
622 }
623 }
624 }
625 } /* next j */
626 } /* next q */
627 } /* next p */
628 } /* next l */
629 } /* next k */
630
631 /* read out matrix minimum */
632 for (j = i + turn + 1; j <= n; j++) {
633 if (evaluate_ext(i, j, i, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
634 j_pos_end = MIN2(n + 1, j + max_interaction_length);
635 for (k = i - 1; k > i_pos_begin; k--) {
636 sk = (k > 1) ? S1[k - 1] : -1;
637 for (l = j + 1; l < j_pos_end; l++) {
638 if (evaluate_ext(k, l, k, l, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
639 type2 = md->pair[S[k]][S[l]];
640 sl = (l < n) ? S1[l + 1] : -1;
641 E = c3[j - 1][max_interaction_length - i + k - 1][l - j] +
642 vrna_E_ext_stem(type2, sk, sl, P) +
643 access_s1[i - k + 1][i] +
644 access_s1[l - j + 1][l];
645
646 if (E < Emin) {
647 Emin = E;
648 k_min = k;
649 l_min = l;
650 j_min = j;
651 }
652 }
653 }
654 }
655 }
656 }
657
658 /*
659 * Only consider hits where the interaction is more
660 * stable than the PK penalty, i.e. the net free
661 * energy is negative
662 */
663 if (Emin < -penalty) {
664 struc = backtrack_XS(fc, k_min, l_min, i, j_min, max_interaction_length, c3);
665
666 dGx = access_s1[i - k_min + 1][i];
667 dGy = access_s1[l_min - j_min + 1][l_min];
668 inter = Emin - dGx - dGy;
669
670 entry = (vrna_pk_plex_t *)vrna_alloc(sizeof(vrna_pk_plex_t));
671 entry->start_5 = k_min;
672 entry->end_5 = i;
673 entry->start_3 = j_min;
674 entry->end_3 = l_min;
675 entry->energy = (double)Emin * 0.01;
676 entry->dG1 = (double)dGx * 0.01;
677 entry->dG2 = (double)dGy * 0.01;
678 entry->dGint = (double)inter * 0.01;
679 entry->structure = struc;
680
681 vrna_heap_insert(storage, entry);
682 }
683 }
684 }
685
686 free_array(c3, n, max_interaction_length);
687
688 return storage;
689 }
690
691
692 PRIVATE char *
backtrack_XS(vrna_fold_compound_t * fc,int k,int l,const int i,const int j,const int max_interaction_length,int *** c3)693 backtrack_XS(vrna_fold_compound_t *fc,
694 int k,
695 int l,
696 const int i,
697 const int j,
698 const int max_interaction_length,
699 int ***c3)
700 {
701 /* backtrack structure going backwards from i, and forwards from j
702 * return structure in bracket notation with & as separator */
703 short *S, *S1;
704 int p, q, type, type2, E, traced, i0, j0, *rtype;
705 char *st1, *st2, *struc;
706 vrna_param_t *P;
707 vrna_md_t *md;
708
709 S = fc->sequence_encoding2;
710 S1 = fc->sequence_encoding;
711 P = fc->params;
712 md = &(P->model_details);
713 rtype = &(md->rtype[0]);
714 st1 = (char *)vrna_alloc(sizeof(char) * (i - k + 2));
715 st1[i - k + 1] = '\0';
716 st2 = (char *)vrna_alloc(sizeof(char) * (l - j + 2));
717 st2[l - j + 1] = '\0';
718
719 i0 = k;
720 j0 = l;
721 while (k <= i && l >= j) {
722 E = c3[j - 1][max_interaction_length - i + k - 1][l - j];
723 traced = 0;
724 st1[k - i0] = '(';
725 st2[l - j] = ')';
726
727 type = md->pair[S[k]][S[l]];
728
729 if (!type)
730 vrna_message_error("backtrack failed in fold duplex bli");
731
732 for (p = k + 1; p <= i; p++) {
733 for (q = l - 1; q >= j; q--) {
734 int LE;
735 if (p - k + l - q - 2 > MAXLOOP)
736 break;
737
738 type2 = md->pair[S[q]][S[p]];
739
740 if (!type2)
741 continue;
742
743 LE = E_IntLoop(p - k - 1,
744 l - q - 1,
745 type,
746 type2,
747 S1[k + 1],
748 S1[l - 1],
749 S1[p - 1],
750 S1[q + 1],
751 P);
752 if (E == c3[j - 1][max_interaction_length - i + p - 1][q - j] + LE) {
753 traced = 1;
754 k = p;
755 l = q;
756 break;
757 }
758 }
759 if (traced)
760 break;
761 }
762 if (!traced) {
763 E -= vrna_E_ext_stem(rtype[type],
764 S1[l - 1],
765 S1[k + 1],
766 P);
767
768 if (E != 0)
769 vrna_message_error("backtrack failed in fold duplex bal");
770 else
771 break;
772 }
773 }
774 struc = (char *)vrna_alloc(sizeof(char) * (k - i0 + 1 + j0 - l + 1 + 2));
775
776 for (p = 0; p <= i - i0; p++)
777 if (!st1[p])
778 st1[p] = '.';
779
780 for (p = 0; p <= j0 - j; p++)
781 if (!st2[p])
782 st2[p] = '.';
783
784 strcpy(struc, st1);
785 strcat(struc, "&");
786 strcat(struc, st2);
787 free(st1);
788 free(st2);
789 return struc;
790 }
791
792
793 PRIVATE
794 int
PKplex_heap_cmp(const void * a,const void * b,void * data)795 PKplex_heap_cmp(const void *a,
796 const void *b,
797 void *data)
798 {
799 vrna_pk_plex_t *p1 = (vrna_pk_plex_t *)a;
800 vrna_pk_plex_t *p2 = (vrna_pk_plex_t *)b;
801
802 if (p1->energy > p2->energy)
803 return 1;
804 else if (p1->energy < p2->energy)
805 return -1;
806
807 return 0;
808 }
809
810
811 /*
812 * ###########################################
813 * # deprecated functions below #
814 *###########################################
815 */
816
817 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
818
819 PUBLIC dupVar *
PKLduplexfold_XS(const char * s1,const int ** access_s1,int penalty,int max_interaction_length,int delta)820 PKLduplexfold_XS(const char *s1,
821 const int **access_s1,
822 int penalty,
823 int max_interaction_length,
824 int delta)
825 {
826 vrna_fold_compound_t *fc;
827 dupVar *hits;
828 default_data scoring_dat;
829 size_t cnt;
830 vrna_heap_t interactions;
831 vrna_pk_plex_t *entry;
832
833 hits = NULL;
834
835 if ((s1) &&
836 (access_s1)) {
837 fc = vrna_fold_compound(s1, NULL, VRNA_OPTION_DEFAULT);
838
839 prepare_PKplex(fc);
840
841 scoring_dat.penalty = -penalty;
842
843 interactions = duplexfold_XS(fc,
844 access_s1,
845 max_interaction_length,
846 &default_pk_plex_penalty,
847 (void *)&scoring_dat);
848
849 cnt = 0;
850 hits = (dupVar *)vrna_alloc(sizeof(dupVar) * (vrna_heap_size(interactions) + 2));
851 while ((entry = vrna_heap_pop(interactions))) {
852 hits[cnt].structure = entry->structure;
853 hits[cnt].tb = entry->start_5;
854 hits[cnt].te = entry->end_5;
855 hits[cnt].qb = entry->start_3;
856 hits[cnt].qe = entry->end_3;
857 hits[cnt].ddG = entry->energy;
858 hits[cnt].dG1 = entry->dG1;
859 hits[cnt].dG2 = entry->dG2;
860 hits[cnt].energy = entry->dGint;
861 hits[cnt].inactive = 0;
862 hits[cnt].processed = 0;
863 free(entry);
864 cnt++;
865 }
866
867 hits[cnt].inactive = 1;
868 hits[cnt].structure = NULL;
869
870 vrna_heap_free(interactions);
871 vrna_fold_compound_free(fc);
872 }
873
874 return hits;
875 }
876
877
878 #endif
879