1 /*
2 * compute the duplex structure of two RNA strands,
3 * allowing only inter-strand base pairs.
4 * see cofold() for computing hybrid structures without
5 * restriction.
6 *
7 * Ivo Hofacker
8 * Vienna RNA package
9 */
10
11 #ifdef HAVE_CONFIG_H
12 #include "config.h"
13 #endif
14
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <math.h>
18 #include <ctype.h>
19 #include <string.h>
20 #include "ViennaRNA/utils/basic.h"
21 #include "ViennaRNA/utils/strings.h"
22 #include "ViennaRNA/params/default.h"
23 #include "ViennaRNA/fold_vars.h"
24 #include "ViennaRNA/snofold.h"
25 #include "ViennaRNA/pair_mat.h"
26 #include "ViennaRNA/params/basic.h"
27 #include "ViennaRNA/snoop.h"
28 #include "ViennaRNA/plotting/probabilities.h"
29 #include "ViennaRNA/plotting/structures.h"
30 /* #include "ViennaRNA/fold.h" */
31 #include "ViennaRNA/duplex.h"
32 #include "ViennaRNA/loops/all.h"
33
34
35 #define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
36 #define NEW_NINIO 1 /* new asymetry penalty */
37
38
39 PRIVATE void
40 encode_seqs(const char *s1,
41 const char *s2);
42
43
44 PRIVATE short *
45 encode_seq(const char *seq);
46
47
48 PRIVATE void
49 find_max_snoop(const char *s1,
50 const char *s2,
51 const int max,
52 const int alignment_length,
53 const int *position,
54 const int delta,
55 const int distance,
56 const int penalty,
57 const int threshloop,
58 const int threshLE,
59 const int threshRE,
60 const int threshDE,
61 const int threshTE,
62 const int threshSE,
63 const int threshD,
64 const int half_stem,
65 const int max_half_stem,
66 const int min_s2,
67 const int max_s2,
68 const int min_s1,
69 const int max_s1,
70 const int min_d1,
71 const int min_d2,
72 const char *name,
73 const int fullStemEnergy);
74
75
76 PRIVATE void
77 find_max_snoop_XS(const char *s1,
78 const char *s2,
79 const int **access_s1,
80 const int max,
81 const int alignment_length,
82 const int *position,
83 const int *position_j,
84 const int delta,
85 const int distance,
86 const int penalty,
87 const int threshloop,
88 const int threshLE,
89 const int threshRE,
90 const int threshDE,
91 const int threshTE,
92 const int threshSE,
93 const int threshD,
94 const int half_stem,
95 const int max_half_stem,
96 const int min_s2,
97 const int max_s2,
98 const int min_s1,
99 const int max_s1,
100 const int min_d1,
101 const int min_d2,
102 const char *name,
103 const int fullStemEnergy);
104
105
106 PRIVATE char *
107 alisnoop_backtrack(int i,
108 int j,
109 const char **s2,
110 int *Duplex_El,
111 int *Duplex_Er,
112 int *Loop_E,
113 int *Loop_D,
114 int *u,
115 int *pscd,
116 int *psct,
117 int *pscg,
118 const int penalty,
119 const int threshloop,
120 const int threshLE,
121 const int threshRE,
122 const int threshDE,
123 const int threshD,
124 const int half_stem,
125 const int max_half_stem,
126 const int min_s2,
127 const int max_s2,
128 const int min_s1,
129 const int max_s1,
130 const int min_d1,
131 const int min_d2,
132 const short **S1,
133 const short **S2);
134
135
136 PRIVATE char *
137 snoop_backtrack(int i,
138 int j,
139 const char *s2,
140 int *Duplex_El,
141 int *Duplex_Er,
142 int *Loop_E,
143 int *Loop_D,
144 int *u,
145 const int penalty,
146 const int threshloop,
147 const int threshLE,
148 const int threshRE,
149 const int threshDE,
150 const int threshD,
151 const int half_stem,
152 const int max_half_stem,
153 const int min_s2,
154 const int max_s2,
155 const int min_s1,
156 const int max_s1,
157 const int min_d1,
158 const int min_d2);
159
160
161 PRIVATE char *
162 snoop_backtrack_XS(int i,
163 int j,
164 const char *s2,
165 int *Duplex_El,
166 int *Duplex_Er,
167 int *Loop_E,
168 int *Loop_D,
169 int *u,
170 const int penalty,
171 const int threshloop,
172 const int threshLE,
173 const int threshRE,
174 const int threshDE,
175 const int threshD,
176 const int half_stem,
177 const int max_half_stem,
178 const int min_s2,
179 const int max_s2,
180 const int min_s1,
181 const int max_s1,
182 const int min_d1,
183 const int min_d2);
184
185
186 PRIVATE int
187 compare(const void *sub1,
188 const void *sub2);
189
190
191 PRIVATE int
192 covscore(const int *types,
193 int n_seq);
194
195
196 PRIVATE short *
197 aliencode_seq(const char *sequence);
198
199
200 PUBLIC int snoop_subopt_sorted = 0; /* from subopt.c, default 0 */
201
202
203 /*@unused@*/
204
205
206 #define MAXLOOP_L 3
207 #define MIN2(A, B) ((A) < (B) ? (A) : (B))
208 #define MAX2(A, B) ((A) > (B) ? (A) : (B))
209 #define ASS 1
210 PRIVATE vrna_param_t *P = NULL;
211
212 PRIVATE int **c = NULL; /* energy array, given that i-j pair */
213 PRIVATE int **r = NULL;
214 PRIVATE int **lc = NULL; /* energy array, given that i-j pair */
215 PRIVATE int **lr = NULL;
216 PRIVATE int **c_fill = NULL;
217 PRIVATE int **r_fill = NULL;
218 PRIVATE int **lpair = NULL;
219
220
221 PRIVATE short *S1 = NULL, *SS1 = NULL, *S2 = NULL, *SS2 = NULL;
222 PRIVATE short *S1_fill = NULL, *SS1_fill = NULL, *S2_fill = NULL, *SS2_fill = NULL;
223 PRIVATE int n1, n2; /* sequence lengths */
224
225 extern int cut_point;
226
227 PRIVATE int delay_free = 0;
228 /*--------------------------------------------------------------------------*/
229
230 snoopT
alisnoopfold(const char ** s1,const char ** s2,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2)231 alisnoopfold(const char **s1,
232 const char **s2,
233 const int penalty,
234 const int threshloop,
235 const int threshLE,
236 const int threshRE,
237 const int threshDE,
238 const int threshD,
239 const int half_stem,
240 const int max_half_stem,
241 const int min_s2,
242 const int max_s2,
243 const int min_s1,
244 const int max_s1,
245 const int min_d1,
246 const int min_d2)
247 {
248 int s, n_seq;
249 int i, j, E, l1, Emin = INF, i_min = 0, j_min = 0;
250 char *struc;
251 snoopT mfe;
252 int *indx;
253 int *mLoop;
254 int *cLoop;
255 folden **foldlist;
256 folden **foldlist_XS;
257 int Duplex_El, Duplex_Er, pscd, psct, pscg;
258 int Loop_D;
259 int u;
260 int Loop_E;
261 short **Sali1, **Sali2;
262 int *type, *type2, *type3;
263 vrna_md_t md;
264
265 Duplex_El = 0;
266 Duplex_Er = 0;
267 Loop_E = 0;
268 Loop_D = 0;
269 pscd = 0;
270 psct = 0;
271 pscg = 0;
272 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
273 n1 = (int)strlen(s1[0]);
274 n2 = (int)strlen(s2[0]);
275
276 for (s = 0; s1[s] != NULL; s++) ;
277 n_seq = s;
278 for (s = 0; s2[s] != NULL; s++) ;
279 if (n_seq != s)
280 vrna_message_error("unequal number of sequences in aliduplexfold()\n");
281
282 set_model_details(&md);
283 if ((!P) || (fabs(P->temperature - temperature) > 1e-6)) {
284 snoupdate_fold_params();
285 if (P)
286 free(P);
287
288 P = vrna_params(&md);
289 make_pair_matrix();
290 }
291
292 c = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
293 r = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
294 for (i = 0; i <= n1; i++) {
295 c[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
296 r[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
297 for (j = n2; j > -1; j--) {
298 c[i][j] = INF;
299 r[i][j] = INF;
300 }
301 }
302 Sali1 = (short **)vrna_alloc((n_seq + 1) * sizeof(short *));
303 Sali2 = (short **)vrna_alloc((n_seq + 1) * sizeof(short *));
304 for (s = 0; s < n_seq; s++) {
305 if ((int)strlen(s1[s]) != n1)
306 vrna_message_error("uneqal seqence lengths");
307
308 if ((int)strlen(s2[s]) != n2)
309 vrna_message_error("uneqal seqence lengths");
310
311 Sali1[s] = aliencode_seq(s1[s]);
312 Sali2[s] = aliencode_seq(s2[s]);
313 }
314 type = (int *)vrna_alloc(n_seq * sizeof(int));
315 type2 = (int *)vrna_alloc(n_seq * sizeof(int));
316 type3 = (int *)vrna_alloc(n_seq * sizeof(int));
317 /* encode_seqs(s1, s2); */
318 for (i = 6; i <= n1 - 5; i++) {
319 int U;
320 U = 0;
321 for (s = 0; s < n_seq; s++)
322 U += Sali1[s][i - 2];
323 U = (U == (n_seq) * 4 ? 1 : 0);
324 for (j = n2 - min_d2; j > min_d1; j--) {
325 int type4, k, l, psc, psc2, psc3;
326 for (s = 0; s < n_seq; s++)
327 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
328 psc = covscore(type, n_seq);
329 for (s = 0; s < n_seq; s++)
330 if (type[s] == 0)
331 type[s] = 7;
332
333 c[i][j] = (psc >= MINPSCORE) ? (n_seq * P->DuplexInit) : INF;
334 if (psc < MINPSCORE)
335 continue;
336
337 if (/* pair[Sali1[i+1]][Sali2[j-1]] && */
338 U && j < max_s1 && j > min_s1 &&
339 j > n2 - max_s2 - max_half_stem &&
340 j < n2 - min_s2 - half_stem) {
341 /*constraint on s2 and i*/
342 folden *temp;
343 temp = foldlist[j + 1];
344 while (temp->next) {
345 int k = temp->k;
346 for (s = 0; s < n_seq; s++) {
347 type2[s] = pair[Sali1[s][i - 3]][Sali2[s][k + 1]];
348 type3[s] = pair[Sali1[s][i - 4]][Sali2[s][k + 1]];
349 }
350 psc2 = covscore(type2, n_seq);
351 psc3 = covscore(type3, n_seq);
352 if (psc2 > MINPSCORE)
353 r[i][j] = MIN2(r[i][j], c[i - 3][k + 1] + temp->energy);
354
355 if (psc3 > MINPSCORE)
356 r[i][j] = MIN2(r[i][j], c[i - 4][k + 1] + temp->energy);
357
358 temp = temp->next;
359 }
360 }
361
362 /* dangle 5'SIDE relative to the mRNA */
363 for (s = 0; s < n_seq; s++)
364 c[i][j] += vrna_E_ext_stem(type[s], Sali1[s][i - 1], Sali2[s][j + 1], P);
365 for (k = i - 1; k > 0 && (i - k) < MAXLOOP_L; k--) {
366 for (l = j + 1; l <= n2; l++) {
367 if (i - k + l - j > 2 * MAXLOOP_L - 2)
368 break;
369
370 if (abs(i - k - l + j) >= ASS)
371 continue;
372
373 for (E = s = 0; s < n_seq; s++) {
374 type4 = pair[Sali1[s][k]][Sali2[s][l]];
375 if (type4 == 0)
376 type4 = 7;
377
378 E += E_IntLoop(i - k - 1, l - j - 1, type4, rtype[type[s]],
379 Sali1[s][k + 1], Sali2[s][l - 1], Sali1[s][i - 1], Sali2[s][j + 1], P);
380 }
381 c[i][j] = MIN2(c[i][j], c[k][l] + E);
382 r[i][j] = MIN2(r[i][j], r[k][l] + E);
383 }
384 }
385 c[i][j] -= psc;
386 r[i][j] -= psc;
387 E = r[i][j];
388 for (s = 0; s < n_seq; s++)
389 E += vrna_E_ext_stem(rtype[type[s]], Sali2[s][j - 1], Sali1[s][i + 1], P);
390 /**
391 *** if (i<n1) E += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
392 *** if (j>1) E += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
393 *** if (type[s]>2) E += P->TerminalAU;
394 **/
395 if (E < Emin) {
396 Emin = E;
397 i_min = i;
398 j_min = j;
399 }
400 }
401 }
402 if (Emin > 0) {
403 printf("no target found under the constraints chosen\n");
404 for (i = 0; i <= n1; i++) {
405 free(r[i]);
406 free(c[i]);
407 }
408 free(c);
409 free(r);
410 for (s = 0; s < n_seq; s++) {
411 free(Sali1[s]);
412 free(Sali2[s]);
413 }
414 free(Sali1);
415 free(Sali2);
416 free(S2);
417 free(SS1);
418 free(SS2);
419 free(type);
420 free(type2);
421 free(type3);
422 mfe.energy = INF;
423 mfe.structure = NULL;
424 return mfe;
425 }
426
427 struc = alisnoop_backtrack(i_min, j_min, (const char **)s2,
428 &Duplex_El, &Duplex_Er, &Loop_E,
429 &Loop_D, &u, &pscd, &psct, &pscg,
430 penalty, threshloop, threshLE,
431 threshRE, threshDE, threshD,
432 half_stem, max_half_stem, min_s2,
433 max_s2, min_s1, max_s1, min_d1,
434 min_d2, (const short **)Sali1, (const short **)Sali2);
435 /* if (i_min<n1-5) i_min++; */
436 /* if (j_min>6 ) j_min--; */
437 l1 = strchr(struc, '&') - struc;
438 mfe.i = i_min - 5;
439 mfe.j = j_min - 5;
440 mfe.u = u - 5;
441 mfe.Duplex_Er = (float)Duplex_Er / 100;
442 mfe.Duplex_El = (float)Duplex_El / 100;
443 mfe.Loop_D = (float)Loop_D / 100;
444 mfe.Loop_E = (float)Loop_E / 100;
445 mfe.energy = (float)Emin / 100;
446 /* mfe.fullStemEnergy = (float) fullStemEnergy/100; */
447 mfe.pscd = pscd;
448 mfe.psct = psct;
449 mfe.structure = struc;
450 for (s = 0; s < n_seq; s++) {
451 free(Sali1[s]);
452 free(Sali2[s]);
453 }
454 free(Sali1);
455 free(Sali2);
456 free(type);
457 free(type2);
458 free(type3);
459
460 if (!delay_free) {
461 for (i = 0; i <= n1; i++) {
462 free(r[i]);
463 free(c[i]);
464 }
465 free(c);
466 free(r);
467 free(S2);
468 free(SS1);
469 free(SS2);
470 }
471
472 return mfe;
473 }
474
475
476 PUBLIC snoopT *
alisnoop_subopt(const char ** s1,const char ** s2,int delta,int w,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int distance,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2)477 alisnoop_subopt(const char **s1,
478 const char **s2,
479 int delta,
480 int w,
481 const int penalty,
482 const int threshloop,
483 const int threshLE,
484 const int threshRE,
485 const int threshDE,
486 const int threshTE,
487 const int threshSE,
488 const int threshD,
489 const int distance,
490 const int half_stem,
491 const int max_half_stem,
492 const int min_s2,
493 const int max_s2,
494 const int min_s1,
495 const int max_s1,
496 const int min_d1,
497 const int min_d2)
498 {
499 short **Sali1, **Sali2;
500 /* printf("%d %d\n", min_s2, max_s2); */
501 int i, j, s, n_seq, n1, n2, E, n_subopt = 0, n_max;
502 char *struc;
503 snoopT mfe;
504 snoopT *subopt;
505 int thresh;
506 int *type;
507 int Duplex_El, Duplex_Er, Loop_E, pscd, psct, pscg;
508 int Loop_D;
509
510 Duplex_El = 0;
511 Duplex_Er = 0;
512 Loop_E = 0;
513 Loop_D = 0;
514 pscd = 0;
515 psct = 0;
516 pscg = 0;
517 int u;
518 u = 0;
519 n_max = 16;
520 subopt = (snoopT *)vrna_alloc(n_max * sizeof(snoopT));
521 delay_free = 1;
522 mfe = alisnoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE, threshD,
523 half_stem, max_half_stem,
524 min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
525 if (mfe.energy > 0) {
526 free(subopt);
527 delay_free = 0;
528 return NULL;
529 }
530
531 thresh = MIN2((int)((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E) * 100 + 0.1 + 410) + delta,
532 threshTE);
533 /* subopt[n_subopt++]=mfe; */
534 free(mfe.structure);
535 n1 = (int)strlen(s1[0]);
536 n2 = (int)strlen(s2[0]);
537 for (s = 0; s1[s] != NULL; s++) ;
538 n_seq = s;
539 Sali1 = (short **)vrna_alloc((n_seq + 1) * sizeof(short *));
540 Sali2 = (short **)vrna_alloc((n_seq + 1) * sizeof(short *));
541 for (s = 0; s < n_seq; s++) {
542 if ((int)strlen(s1[s]) != n1)
543 vrna_message_error("uneqal seqence lengths");
544
545 if ((int)strlen(s2[s]) != n2)
546 vrna_message_error("uneqal seqence lengths");
547
548 Sali1[s] = aliencode_seq(s1[s]);
549 Sali2[s] = aliencode_seq(s2[s]);
550 }
551 Sali1[n_seq] = NULL;
552 Sali2[n_seq] = NULL;
553 type = (int *)vrna_alloc(n_seq * sizeof(int));
554 for (i = n1; i > 1; i--) {
555 for (j = 1; j <= n2; j++) {
556 int ii, jj, Ed, psc, skip;
557 for (s = 0; s < n_seq; s++)
558 type[s] = pair[Sali2[s][j]][Sali1[s][i]];
559 psc = covscore(type, n_seq);
560 for (s = 0; s < n_seq; s++)
561 if (type[s] == 0)
562 type[s] = 7;
563
564 if (psc < MINPSCORE)
565 continue;
566
567 E = Ed = r[i][j];
568 for (s = 0; s < n_seq; s++) {
569 /* if (i<n1-5) Ed += P->dangle3[type[s]][Sali1[s][i+1]]; */
570 /* if (j>6) Ed += P->dangle5[type[s]][Sali2[s][j-1]]; */
571 if (type[s] > 2)
572 Ed += P->TerminalAU;
573 }
574 if (Ed > thresh)
575 continue;
576
577 /* too keep output small, remove hits that are dominated by a
578 * better one close (w) by. For simplicity we do test without
579 * adding dangles, which is slightly inaccurate.
580 */
581 w = 1;
582 for (skip = 0, ii = MAX2(i - w, 1); (ii <= MIN2(i + w, n1)) && type; ii++) {
583 for (jj = MAX2(j - w, 1); jj <= MIN2(j + w, n2); jj++)
584 if (r[ii][jj] < E) {
585 skip = 1;
586 break;
587 }
588 }
589 if (skip)
590 continue;
591
592 psct = 0;
593 pscg = 0;
594 struc = alisnoop_backtrack(i,
595 j,
596 s2,
597 &Duplex_El,
598 &Duplex_Er,
599 &Loop_E,
600 &Loop_D,
601 &u,
602 &pscd,
603 &psct,
604 &pscg,
605 penalty,
606 threshloop,
607 threshLE,
608 threshRE,
609 threshDE,
610 threshD,
611 half_stem,
612 max_half_stem,
613 min_s2,
614 max_s2,
615 min_s1,
616 max_s1,
617 min_d1,
618 min_d2,
619 (const short int **)Sali1,
620 (const int short **)Sali2);
621
622 if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
623 (Duplex_Er + Duplex_El) > threshDE ||
624 (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
625 (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
626 /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
627 /* " Duplex_Er + Duplex_El %d threshDE %d \n" */
628 /* " Duplex_Er + Duplex_El + Loop_E %d threshTE %d \n" */
629 /* " Duplex_Er + Duplex_El + Loop_E + Loop_D %d threshSE %d \n", */
630 /* Duplex_Er , threshRE , Duplex_El ,threshLE, */
631 /* Duplex_Er + Duplex_El, threshDE, */
632 /* Duplex_Er + Duplex_El+ Loop_E , threshTE, */
633 /* Duplex_Er + Duplex_El+ Loop_E + Loop_D, threshSE); */
634 Duplex_Er = 0;
635 Duplex_El = 0;
636 Loop_E = 0;
637 Loop_D = 0;
638 u = 0,
639 free(struc);
640 continue;
641 }
642
643 if (n_subopt + 1 >= n_max) {
644 n_max *= 2;
645 subopt = (snoopT *)vrna_realloc(subopt, n_max * sizeof(snoopT));
646 }
647
648 subopt[n_subopt].i = i - 5;
649 subopt[n_subopt].j = j - 5;
650 subopt[n_subopt].u = u - 5;
651 subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
652 subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
653 subopt[n_subopt].Loop_E = Loop_E * 0.01;
654 subopt[n_subopt].Loop_D = Loop_D * 0.01;
655 subopt[n_subopt].energy = (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) * 0.01;
656 subopt[n_subopt].pscd = pscd * 0.01;
657 subopt[n_subopt].psct = -psct * 0.01;
658 subopt[n_subopt++].structure = struc;
659
660 /* i=u; */
661 Duplex_Er = 0;
662 Duplex_El = 0;
663 Loop_E = 0;
664 Loop_D = 0;
665 u = 0;
666 pscd = 0;
667 psct = 0;
668 }
669 }
670
671 for (i = 0; i <= n1; i++) {
672 free(c[i]);
673 free(r[i]);
674 }
675 free(c);
676 free(r);
677 for (s = 0; s < n_seq; s++) {
678 free(Sali1[s]);
679 free(Sali2[s]);
680 }
681 free(Sali1);
682 free(Sali2);
683 free(type);
684
685 if (snoop_subopt_sorted)
686 qsort(subopt, n_subopt, sizeof(snoopT), compare);
687
688 subopt[n_subopt].i = 0;
689 subopt[n_subopt].j = 0;
690 subopt[n_subopt].structure = NULL;
691 return subopt;
692 }
693
694
695 PRIVATE char *
alisnoop_backtrack(int i,int j,const char ** snoseq,int * Duplex_El,int * Duplex_Er,int * Loop_E,int * Loop_D,int * u,int * pscd,int * psct,int * pscg,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const short ** Sali1,const short ** Sali2)696 alisnoop_backtrack(int i,
697 int j,
698 const char **snoseq,
699 int *Duplex_El,
700 int *Duplex_Er,
701 int *Loop_E,
702 int *Loop_D,
703 int *u,
704 int *pscd,
705 int *psct,
706 int *pscg,
707 const int penalty,
708 const int threshloop,
709 const int threshLE,
710 const int threshRE,
711 const int threshDE,
712 const int threshD,
713 const int half_stem,
714 const int max_half_stem,
715 const int min_s2,
716 const int max_s2,
717 const int min_s1,
718 const int max_s1,
719 const int min_d1,
720 const int min_d2,
721 const short **Sali1,
722 const short **Sali2)
723 {
724 /* backtrack structure going backwards from i, and forwards from j
725 * return structure in bracket notation with & as separator */
726 int k, l, *type, *type2, *type3, type4, E, traced, i0, j0, s, n_seq, psc;
727 int traced_r = 0; /* flag for following backtrack in c or r */
728 char *st1, *st2, *struc;
729 char *struc_loop;
730
731 n1 = (int)Sali1[0][0];
732 n2 = (int)Sali2[0][0];
733
734 for (s = 0; Sali1[s] != NULL; s++) ;
735 n_seq = s;
736 for (s = 0; Sali2[s] != NULL; s++) ;
737 if (n_seq != s)
738 vrna_message_error("unequal number of sequences in alibacktrack()\n");
739
740 st1 = (char *)vrna_alloc(sizeof(char) * (n1 + 1));
741 st2 = (char *)vrna_alloc(sizeof(char) * (n2 + 1));
742 type = (int *)vrna_alloc(n_seq * sizeof(int));
743 type2 = (int *)vrna_alloc(n_seq * sizeof(int));
744 type3 = (int *)vrna_alloc(n_seq * sizeof(int));
745 int *indx;
746 int *mLoop;
747 int *cLoop;
748 folden **foldlist, **foldlist_XS;
749 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
750 i0 = i;
751 j0 = j; /* MIN2(i+1,n1); j0=MAX2(j-1,1);!modified */
752 for (s = 0; s < n_seq; s++) {
753 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
754 if (type[s] == 0)
755 type[s] = 7;
756
757 *Duplex_Er +=
758 vrna_E_ext_stem(rtype[type[s]], (j > 1) ? Sali2[s][j - 1] : -1, (i < n1) ? Sali1[s][i + 1] : -1, P);
759 /**
760 *** if (i<n1) *Duplex_Er += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
761 *** if (j>1) *Duplex_Er += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
762 *** if (type[s]>2) *Duplex_Er += P->TerminalAU;
763 **/
764 }
765 while (i > 0 && j <= n2 - min_d2) {
766 if (!traced_r) {
767 E = r[i][j];
768 traced = 0;
769 st1[i - 1] = '<';
770 st2[j - 1] = '>';
771 for (s = 0; s < n_seq; s++)
772 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
773 psc = covscore(type, n_seq);
774 for (s = 0; s < n_seq; s++)
775 if (type[s] == 0)
776 type[s] = 7;
777
778 E += psc;
779 *pscd += psc;
780 for (k = i - 1; k > 0 && (i - k) < MAXLOOP_L; k--) {
781 for (l = j + 1; l <= n2; l++) {
782 int LE;
783 if (i - k + l - j > 2 * MAXLOOP_L - 2)
784 break;
785
786 if (abs(i - k - l + j) >= ASS)
787 continue;
788
789 for (s = LE = 0; s < n_seq; s++) {
790 type4 = pair[Sali1[s][k]][Sali2[s][l]];
791 if (type4 == 0)
792 type4 = 7;
793
794 LE +=
795 E_IntLoop(i - k - 1,
796 l - j - 1,
797 type4,
798 rtype[type[s]],
799 Sali1[s][k + 1],
800 Sali2[s][l - 1],
801 Sali1[s][i - 1],
802 Sali2[s][j + 1],
803 P);
804 }
805 if (E == r[k][l] + LE) {
806 traced = 1;
807 i = k;
808 j = l;
809 *Duplex_Er += LE;
810 break;
811 }
812 }
813 if (traced)
814 break;
815 }
816 if (!traced) {
817 int U = 0;
818 for (s = 0; s < n_seq; s++)
819 U += Sali1[s][i - 2];
820 U = (U == (n_seq) * 4 ? 1 : 0);
821 if (/* pair[Sali1[i+1]][Sali2[j-1]] && */ /* only U's are allowed */
822 U && j < max_s1 && j > min_s1 &&
823 j > n2 - max_s2 - max_half_stem &&
824 j < n2 - min_s2 - half_stem) {
825 int min_k, max_k;
826 max_k = MIN2(n2 - min_s2, j + max_half_stem + 1);
827 min_k = MAX2(j + half_stem + 1, n2 - max_s2);
828 folden *temp;
829 temp = foldlist[j + 1];
830 while (temp->next) {
831 int psc2, psc3;
832 int k = temp->k;
833 for (s = 0; s < n_seq; s++) {
834 type2[s] = pair[Sali1[s][i - 3]][Sali2[s][k + 1]];
835 type3[s] = pair[Sali1[s][i - 4]][Sali2[s][k + 1]];
836 }
837 psc2 = covscore(type2, n_seq);
838 psc3 = covscore(type3, n_seq);
839 if (psc2 > MINPSCORE /*&& pair[Sali1[i-4]][Sali2[k+2]]*/) {
840 /* introduce structure from RNAfold */
841 if (E == c[i - 3][k + 1] + temp->energy) {
842 *Loop_E = temp->energy;
843 st1[i - 3] = '|';
844 *u = i - 2;
845 int a, b;
846 /* int fix_ij=indx[k-1+1]+j+1; */
847 for (a = 0; a < MISMATCH; a++) {
848 for (b = 0; b < MISMATCH; b++) {
849 int ij = indx[k - 1 - a + 1] + j + 1 + b;
850 if (cLoop[ij] == temp->energy) {
851 /* int bla; */
852 struc_loop = alisnobacktrack_fold_from_pair(snoseq,
853 j + 1 + b,
854 k - a - 1 + 1,
855 psct);
856 a = INF;
857 b = INF;
858 }
859 }
860 }
861 traced = 1;
862 traced_r = 1;
863 i = i - 3;
864 j = k + 1;
865 break;
866 }
867 }
868
869 if (psc3 > MINPSCORE /*&& pair[Sali1[i-5]][Sali2[k+2]]*/) {
870 /* introduce structure from RNAfold */
871 if (E == c[i - 4][k + 1] + temp->energy) {
872 *Loop_E = temp->energy;
873 st1[i - 3] = '|';
874 *u = i - 2;
875 int a, b;
876 /* int fix_ij=indx[k-1+1]+j+1; */
877 for (a = 0; a < MISMATCH; a++) {
878 for (b = 0; b < MISMATCH; b++) {
879 int ij = indx[k - 1 - a + 1] + j + 1 + b;
880 if (cLoop[ij] == temp->energy) {
881 /* int bla; */
882 struc_loop = alisnobacktrack_fold_from_pair(snoseq,
883 j + 1 + b,
884 k - a - 1 + 1,
885 psct);
886 a = INF;
887 b = INF;
888 }
889 }
890 }
891 traced = 1;
892 traced_r = 1;
893 i = i - 4;
894 j = k + 1;
895 break;
896 }
897 } /* else if */
898
899 temp = temp->next;
900 } /* while temp-> next */
901 } /* test on j */
902 } /* traced? */
903 } /* traced_r? */
904 else {
905 E = c[i][j];
906 traced = 0;
907 st1[i - 1] = '<';
908 st2[j - 1] = '>';
909 for (s = 0; s < n_seq; s++)
910 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
911 psc = covscore(type, n_seq);
912 for (s = 0; s < n_seq; s++)
913 if (type[s] == 0)
914 type[s] = 7;
915
916 E += psc;
917 *pscd += psc;
918 if (!type)
919 vrna_message_error("backtrack failed in fold duplex c");
920
921 for (k = i - 1; (i - k) < MAXLOOP_L; k--) {
922 for (l = j + 1; l <= n2; l++) {
923 int LE;
924 if (i - k + l - j > 2 * MAXLOOP_L - 2)
925 break;
926
927 if (abs(i - k - l + j) >= ASS)
928 continue;
929
930 for (s = LE = 0; s < n_seq; s++) {
931 type4 = pair[Sali1[s][k]][Sali2[s][l]];
932 if (type4 == 0)
933 type4 = 7;
934
935 LE +=
936 E_IntLoop(i - k - 1,
937 l - j - 1,
938 type4,
939 rtype[type[s]],
940 Sali1[s][k + 1],
941 Sali2[s][l - 1],
942 Sali1[s][i - 1],
943 Sali2[s][j + 1],
944 P);
945 }
946 if (E == c[k][l] + LE) {
947 traced = 1;
948 i = k;
949 j = l;
950 *Duplex_El += LE;
951 break;
952 }
953 }
954 if (traced)
955 break;
956 }
957 }
958
959 if (!traced) {
960 for (s = 0; s < n_seq; s++) {
961 int correction;
962 correction = vrna_E_ext_stem(type[s],
963 (i > 1) ? Sali1[s][i - 1] : -1,
964 (j < n2) ? Sali2[s][j + 1] : -1,
965 P);
966 *Duplex_El += correction;
967 E -= correction;
968 /**
969 *** if (i>1) {E -= P->dangle5[type[s]][Sali1[s][i-1]]; *Duplex_El +=P->dangle5[type[s]][Sali1[s][i-1]];}
970 *** if (j<n2) {E -= P->dangle3[type[s]][Sali2[s][j+1]]; *Duplex_El +=P->dangle3[type[s]][Sali2[s][j+1]];}
971 *** if (type[s]>2) {E -= P->TerminalAU; *Duplex_El +=P->TerminalAU;}
972 **/
973 }
974 if (E != n_seq * P->DuplexInit)
975 vrna_message_error("backtrack failed in fold duplex end");
976 else
977 break;
978 }
979 }
980 /* if (i>1) i--; */
981 /* if (j<n2) j++; */
982 /* struc = (char *) vrna_alloc(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
983 struc = (char *)vrna_alloc(i0 - i + 1 + n2 - 1 + 1 + 2); /* declare final duplex structure */
984 char *struc2;
985 struc2 = (char *)vrna_alloc(n2 + 1);
986 /* char * struct_const; */
987 for (k = MAX2(i, 1); k <= i0; k++)
988 if (!st1[k - 1])
989 st1[k - 1] = '.';
990
991 /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
992 /* char * struct_const; */
993 /* struct_const = (char *) vrna_alloc(sizeof(char)*(n2+1)); */
994 for (k = 1; k <= n2; k++) {
995 if (!st2[k - 1])
996 st2[k - 1] = struc_loop[k - 1]; /* '.'; */
997
998 struc2[k - 1] = st2[k - 1]; /* '.'; */
999 /* if (k>=j0 && k<=j){ */
1000 /* struct_const[k-1]='x'; */
1001 /* } */
1002 /* else{ */
1003 /* if(k<j0) {struct_const[k-1]='<';} */
1004 /* if(k>j) {struct_const[k-1]='>';} */
1005 /* } */
1006 }
1007
1008 /* char duplexseq_1[j0+1]; */
1009 /* char duplexseq_2[n2-j+3]; */
1010 if (j < n2) {
1011 char **duplexseq_1, **duplexseq_2;
1012 duplexseq_1 = (char **)vrna_alloc((n_seq + 1) * sizeof(char *));
1013 duplexseq_2 = (char **)vrna_alloc((n_seq + 1) * sizeof(char *));
1014 for (s = 0; s < n_seq; s++) {
1015 duplexseq_1[s] = (char *)vrna_alloc((j0) * sizeof(char)); /* modfied j0+1 */
1016 duplexseq_2[s] = (char *)vrna_alloc((n2 - j + 2) * sizeof(char)); /* modified j+3 */
1017 strncpy(duplexseq_1[s], snoseq[s], j0 - 1); /* modified j0 */
1018 strcpy(duplexseq_2[s], snoseq[s] + j); /* modified j-1 */
1019 duplexseq_1[s][j0 - 1] = '\0'; /* modified j0 */
1020 duplexseq_2[s][n2 - j + 1] = '\0'; /* modified j+2 */
1021 }
1022 duplexseq_1[n_seq] = NULL;
1023 duplexseq_2[n_seq] = NULL;
1024 duplexT temp;
1025 temp = aliduplexfold((const char **)duplexseq_1, (const char **)duplexseq_2);
1026 *Loop_D = MIN2(0, -410 + (int)100 * temp.energy * n_seq);
1027 if (*Loop_D) {
1028 int l1, ibegin, iend, jbegin, jend;
1029 l1 = strchr(temp.structure, '&') - temp.structure;
1030 ibegin = temp.i - l1;
1031 iend = temp.i - 1;
1032 jbegin = temp.j;
1033 jend = temp.j + (int)strlen(temp.structure) - l1 - 2 - 1;
1034 for (k = ibegin + 1; k <= iend + 1; k++)
1035 struc2[k - 1] = temp.structure[k - ibegin - 1];
1036 for (k = jbegin + j; k <= jend + j; k++)
1037 struc2[k - 1] = temp.structure[l1 + k - j - jbegin + 1];
1038 }
1039
1040 for (s = 0; s < n_seq; s++) {
1041 free(duplexseq_1[s]);
1042 free(duplexseq_2[s]);
1043 }
1044 free(duplexseq_1);
1045 free(duplexseq_2);
1046 free(temp.structure);
1047 }
1048
1049 strcpy(struc, st1 + MAX2(i - 1, 0));
1050 strcat(struc, "&");
1051 /* strcat(struc, st2); */
1052 strncat(struc, struc2 + 5, (int)strlen(struc2) - 10);
1053 free(struc2);
1054 free(struc_loop);
1055 free(st1);
1056 free(st2);
1057 free(type);
1058 free(type2);
1059 free(type3);
1060 /* free_arrays(); */
1061 return struc;
1062 }
1063
1064
1065 void
Lsnoop_subopt(const char * s1,const char * s2,int delta,int w,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int distance,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const int alignment_length,const char * name,const int fullStemEnergy)1066 Lsnoop_subopt(const char *s1,
1067 const char *s2,
1068 int delta,
1069 int w,
1070 const int penalty,
1071 const int threshloop,
1072 const int threshLE,
1073 const int threshRE,
1074 const int threshDE,
1075 const int threshTE,
1076 const int threshSE,
1077 const int threshD,
1078 const int distance,
1079 const int half_stem,
1080 const int max_half_stem,
1081 const int min_s2,
1082 const int max_s2,
1083 const int min_s1,
1084 const int max_s1,
1085 const int min_d1,
1086 const int min_d2,
1087 const int alignment_length,
1088 const char *name,
1089 const int fullStemEnergy)
1090 {
1091 int min_colonne = INF;
1092 int max_pos;
1093 int max;
1094
1095 max = INF;
1096
1097 /* int temp; */
1098 /* int nsubopt=10; */
1099 n1 = (int)strlen(s1);
1100 n2 = (int)strlen(s2);
1101 int *position;
1102 position = (int *)vrna_alloc((n1 + 3) * sizeof(int));
1103
1104
1105 /* int Eminj, Emin_l; */
1106 int i, j; /* l1, Emin=INF, i_min=0, j_min=0; */
1107 /* char *struc; */
1108 /* snoopT mfe; */
1109 int *indx;
1110 int *mLoop;
1111 int *cLoop;
1112 folden **foldlist, **foldlist_XS;
1113 int Duplex_El, Duplex_Er;
1114 int Loop_D;
1115 /* int u; */
1116 int Loop_E;
1117 vrna_md_t md;
1118
1119 Duplex_El = 0;
1120 Duplex_Er = 0;
1121 Loop_E = 0, Loop_D = 0;
1122 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
1123 set_model_details(&md);
1124 if ((!P) || (fabs(P->temperature - temperature) > 1e-6)) {
1125 snoupdate_fold_params();
1126 if (P)
1127 free(P);
1128
1129 P = vrna_params(&md);
1130 make_pair_matrix();
1131 }
1132
1133 lc = (int **)vrna_alloc(sizeof(int *) * (5));
1134 lr = (int **)vrna_alloc(sizeof(int *) * (5));
1135 for (i = 0; i < 5; i++) {
1136 lc[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1137 lr[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1138 for (j = n2; j > -1; j--) {
1139 lc[i][j] = INF;
1140 lr[i][j] = INF;
1141 }
1142 }
1143 encode_seqs(s1, s2);
1144 for (i = 1; i <= n1; i++) {
1145 int idx = i % 5;
1146 int idx_1 = (i - 1) % 5;
1147 int idx_2 = (i - 2) % 5;
1148 int idx_3 = (i - 3) % 5;
1149 int idx_4 = (i - 4) % 5;
1150 for (j = n2 - min_d2; j > min_d1; j--) {
1151 int type, type2, k;
1152 type = pair[S1[i]][S2[j]];
1153 lc[idx][j] = (type) ? P->DuplexInit + 2 * penalty : INF;
1154 lr[idx][j] = INF;
1155 if (!type)
1156 continue;
1157
1158 if ( /*pair[S1[i+1]][S2[j-1]] && check that we have a solid base stack after the mLoop */
1159 j < max_s1 && j > min_s1 &&
1160 j > n2 - max_s2 - max_half_stem &&
1161 j < n2 - min_s2 - half_stem && S1[i - 2] == 4) {
1162 /*constraint on s2 and i*/
1163 int min_k, max_k;
1164 max_k = MIN2(n2 - min_s2, j + max_half_stem + 1);
1165 min_k = MAX2(j + half_stem + 1, n2 - max_s2);
1166 for (k = min_k; k <= max_k; k++) {
1167 if (mLoop[indx[k - 1] + j + 1] < 0) {
1168 }
1169
1170 if (pair[S1[i - 3]][S2[k]] /*genau zwei ungepaarte nucleotiden --NU--*/
1171 && mLoop[indx[k - 1] + j + 1] < threshloop)
1172 lr[idx][j] = MIN2(lr[idx][j], lc[idx_3][k] + mLoop[indx[k - 1] + j + 1]);
1173 else if (pair[S1[i - 4]][S2[k]] && mLoop[indx[k - 1] + j + 1] < threshloop) /*--NUN--*/
1174 lr[idx][j] = MIN2(lr[idx][j], lc[idx_4][k] + mLoop[indx[k - 1] + j + 1]);
1175 }
1176 }
1177
1178 /* dangle 5'SIDE relative to the mRNA */
1179 lc[idx][j] += vrna_E_ext_stem(type, (i > 1) ? SS1[i - 1] : -1, (j < n2) ? SS2[j + 1] : -1, P);
1180 /**
1181 *** if (i>1) lc[idx][j] += P->dangle5[type][SS1[i-1]];
1182 *** if (j<n2) lc[idx][j] += P->dangle3[type][SS2[j+1]];
1183 *** if (type>2) lc[idx][j] += P->TerminalAU;
1184 **/
1185
1186 if (j < n2 && i > 1) {
1187 type2 = pair[S1[i - 1]][S2[j + 1]];
1188 if (type2 > 0) {
1189 lc[idx][j] =
1190 MIN2(lc[idx_1][j + 1] +
1191 E_IntLoop(0, 0, type2, rtype[type], SS1[i], SS2[j], SS1[i - 1], SS2[j + 1],
1192 P) + 2 * penalty,
1193 lc[idx][j]);
1194 lr[idx][j] =
1195 MIN2(lr[idx_1][j + 1] +
1196 E_IntLoop(0, 0, type2, rtype[type], SS1[i], SS2[j], SS1[i - 1], SS2[j + 1],
1197 P) + 2 * penalty,
1198 lr[idx][j]);
1199 }
1200 }
1201
1202 if (j < n2 - 1 && i > 2) {
1203 type2 = pair[S1[i - 2]][S2[j + 2]];
1204 if (type2 > 0) {
1205 lc[idx][j] =
1206 MIN2(lc[idx_2][j + 2] +
1207 E_IntLoop(1, 1, type2, rtype[type], SS1[i - 1], SS2[j + 1], SS1[i - 1], SS2[j + 1],
1208 P) + 4 * penalty,
1209 lc[idx][j]);
1210 lr[idx][j] =
1211 MIN2(lr[idx_2][j + 2] +
1212 E_IntLoop(1, 1, type2, rtype[type], SS1[i - 1], SS2[j + 1], SS1[i - 1], SS2[j + 1],
1213 P) + 4 * penalty,
1214 lr[idx][j]);
1215 }
1216 }
1217
1218 if (j < n2 - 2 && i > 3) {
1219 type2 = pair[S1[i - 3]][S2[j + 3]];
1220 if (type2 > 0) {
1221 lc[idx][j] =
1222 MIN2(lc[idx_3][j + 3] +
1223 E_IntLoop(2, 2, type2, rtype[type], SS1[i - 2], SS2[j + 2], SS1[i - 1], SS2[j + 1],
1224 P) + 6 * penalty,
1225 lc[idx][j]);
1226 lr[idx][j] =
1227 MIN2(lr[idx_3][j + 3] +
1228 E_IntLoop(2, 2, type2, rtype[type], SS1[i - 2], SS2[j + 2], SS1[i - 1], SS2[j + 1],
1229 P) + 6 * penalty,
1230 lr[idx][j]);
1231 }
1232 }
1233
1234 /**
1235 *** (type>2?P->TerminalAU:0)+(i<(n1)?P->dangle3[rtype[type]][SS1[i+1]]+penalty:0)+(j>1?P->dangle5[rtype[type]][SS2[j-1]]+penalty:0)
1236 **/
1237 min_colonne =
1238 MIN2(lr[idx][j] +
1239 vrna_E_ext_stem(rtype[type], (j > 1) ? SS2[j - 1] : -1, (i < n1) ? SS1[i + 1] : -1, P),
1240 min_colonne);
1241 }
1242 position[i] = min_colonne;
1243 if (max >= min_colonne) {
1244 max = min_colonne;
1245 max_pos = i;
1246 }
1247
1248 min_colonne = INF;
1249 }
1250
1251 free(S1);
1252 free(S2);
1253 free(SS1);
1254 free(SS2);
1255 if (max < threshTE) {
1256 find_max_snoop(s1,
1257 s2,
1258 max,
1259 alignment_length,
1260 position,
1261 delta,
1262 distance,
1263 penalty,
1264 threshloop,
1265 threshLE,
1266 threshRE,
1267 threshDE,
1268 threshTE,
1269 threshSE,
1270 threshD,
1271 half_stem,
1272 max_half_stem,
1273 min_s2,
1274 max_s2,
1275 min_s1,
1276 max_s1,
1277 min_d1,
1278 min_d2,
1279 name,
1280 fullStemEnergy);
1281 }
1282
1283 for (i = 1; i < 5; i++) {
1284 free(lc[i]);
1285 free(lr[i]);
1286 }
1287 free(lc[0]);
1288 free(lr[0]);
1289 free(lc);
1290 free(lr);
1291 free(position);
1292 }
1293
1294
1295 void
Lsnoop_subopt_list(const char * s1,const char * s2,int delta,int w,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int distance,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const int alignment_length,const char * name,const int fullStemEnergy)1296 Lsnoop_subopt_list(const char *s1,
1297 const char *s2,
1298 int delta,
1299 int w,
1300 const int penalty,
1301 const int threshloop,
1302 const int threshLE,
1303 const int threshRE,
1304 const int threshDE,
1305 const int threshTE,
1306 const int threshSE,
1307 const int threshD,
1308 const int distance,
1309 const int half_stem,
1310 const int max_half_stem,
1311 const int min_s2,
1312 const int max_s2,
1313 const int min_s1,
1314 const int max_s1,
1315 const int min_d1,
1316 const int min_d2,
1317 const int alignment_length,
1318 const char *name,
1319 const int fullStemEnergy)
1320 {
1321 int min_colonne = INF;
1322 int max_pos;
1323 int max;
1324
1325 max = INF;
1326
1327 /* int temp; */
1328 /* int nsubopt=10; */
1329 n1 = (int)strlen(s1);
1330 n2 = (int)strlen(s2);
1331 int *position;
1332 position = (int *)vrna_alloc((n1 + 3) * sizeof(int));
1333
1334
1335 /* int Eminj, Emin_l; */
1336 int i, j;/* l1, Emin=INF, i_min=0, j_min=0; */
1337 /* char *struc; */
1338 /* snoopT mfe; */
1339 int *indx;
1340 int *mLoop;
1341 int *cLoop;
1342 folden **foldlist, **foldlist_XS;
1343 int Duplex_El, Duplex_Er;
1344 int Loop_D;
1345 /* int u; */
1346 int Loop_E;
1347 vrna_md_t md;
1348
1349 Duplex_El = 0;
1350 Duplex_Er = 0;
1351 Loop_E = 0, Loop_D = 0;
1352 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
1353
1354 set_model_details(&md);
1355 if ((!P) || (fabs(P->temperature - temperature) > 1e-6)) {
1356 snoupdate_fold_params();
1357 if (P)
1358 free(P);
1359
1360 P = vrna_params(&md);
1361 make_pair_matrix();
1362 }
1363
1364 lpair = (int **)vrna_alloc(sizeof(int *) * (6));
1365 lc = (int **)vrna_alloc(sizeof(int *) * (6));
1366 lr = (int **)vrna_alloc(sizeof(int *) * (6));
1367 for (i = 0; i < 6; i++) {
1368 lc[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1369 lr[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1370 lpair[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1371 for (j = n2; j > -1; j--) {
1372 lc[i][j] = INF;
1373 lr[i][j] = INF;
1374 lpair[i][j] = 0;
1375 }
1376 }
1377 encode_seqs(s1, s2);
1378 int lim_maxj = n2 - min_d2;
1379 int lim_minj = min_d1;
1380 int lim_maxi = n1;
1381 for (i = 5; i <= lim_maxi; i++) {
1382 int idx = i % 5;
1383 int idx_1 = (i - 1) % 5;
1384 int idx_2 = (i - 2) % 5;
1385 int idx_3 = (i - 3) % 5;
1386 int idx_4 = (i - 4) % 5;
1387
1388 for (j = lim_maxj; j > lim_minj; j--) {
1389 int type, type2;/* E, k,l; */
1390 type = pair[S1[i]][S2[j]];
1391 lpair[idx][j] = type;
1392 lc[idx][j] = (type) ? P->DuplexInit + 2 * penalty : INF;
1393 lr[idx][j] = INF;
1394 if (!type)
1395 continue;
1396
1397 if ( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
1398 j < max_s1 && j > min_s1 &&
1399 j > n2 - max_s2 - max_half_stem &&
1400 j < n2 - min_s2 - half_stem && S1[i - 2] == 4) {
1401 /*constraint on s2 and i*/
1402 int min_k, max_k;
1403 max_k = MIN2(n2 - min_s2, j + max_half_stem + 1);
1404 min_k = MAX2(j + half_stem + 1, n2 - max_s2);
1405 folden *temp;
1406 temp = foldlist[j + 1];
1407 while (temp->next) {
1408 int k = temp->k;
1409 /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
1410 if (lpair[idx_3][k + 1] /*&& lpair[idx_4][k+2]*/)
1411 lr[idx][j] = MIN2(lr[idx][j], lc[idx_3][k + 1] + temp->energy); /*--NU--*/
1412
1413 /*else*/ if (lpair[idx_4][k + 1]) /*--NUN--*/
1414 lr[idx][j] = MIN2(lr[idx][j], lc[idx_4][k + 1] + temp->energy);
1415
1416 /* } */
1417 temp = temp->next;
1418 }
1419 }
1420
1421 /* dangle 5'SIDE relative to the mRNA */
1422 lc[idx][j] += vrna_E_ext_stem(type, SS1[i - 1], SS2[j + 1], P);
1423 /**
1424 *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
1425 *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
1426 *** if (type>2) lc[idx][j] += P->TerminalAU;
1427 **/
1428 /* if(j<n2 && i>1){ */
1429 /* type2=pair[S1[i-1]][S2[j+1]]; */
1430 type2 = lpair[idx_1][j + 1];
1431 if (type2 > 0) {
1432 lc[idx][j] =
1433 MIN2(lc[idx_1][j + 1] +
1434 E_IntLoop(0, 0, type2, rtype[type], SS1[i], SS2[j], SS1[i - 1], SS2[j + 1],
1435 P) + 2 * penalty,
1436 lc[idx][j]);
1437 lr[idx][j] =
1438 MIN2(lr[idx_1][j + 1] +
1439 E_IntLoop(0, 0, type2, rtype[type], SS1[i], SS2[j], SS1[i - 1], SS2[j + 1],
1440 P) + 2 * penalty,
1441 lr[idx][j]);
1442 }
1443
1444 /* } */
1445 /* if(j<n2-1 && i>2){ */
1446 /* type2=pair[S1[i-2]][S2[j+2]]; */
1447 type2 = lpair[idx_2][j + 2];
1448 if (type2 > 0) {
1449 lc[idx][j] =
1450 MIN2(lc[idx_2][j + 2] +
1451 E_IntLoop(1, 1, type2, rtype[type], SS1[i - 1], SS2[j + 1], SS1[i - 1], SS2[j + 1],
1452 P),
1453 lc[idx][j]);
1454 lr[idx][j] =
1455 MIN2(lr[idx_2][j + 2] +
1456 E_IntLoop(1, 1, type2, rtype[type], SS1[i - 1], SS2[j + 1], SS1[i - 1], SS2[j + 1],
1457 P),
1458 lr[idx][j]);
1459 /* } */
1460 }
1461
1462 /* if(j<n2-2 && i>3){ */
1463 /* type2 = pair[S1[i-3]][S2[j+3]]; */
1464 type2 = lpair[idx_3][j + 3];
1465 if (type2 > 0) {
1466 lc[idx][j] =
1467 MIN2(lc[idx_3][j + 3] +
1468 E_IntLoop(2, 2, type2, rtype[type], SS1[i - 2], SS2[j + 2], SS1[i - 1], SS2[j + 1],
1469 P) + 6 * penalty,
1470 lc[idx][j]);
1471 lr[idx][j] =
1472 MIN2(lr[idx_3][j + 3] +
1473 E_IntLoop(2, 2, type2, rtype[type], SS1[i - 2], SS2[j + 2], SS1[i - 1], SS2[j + 1],
1474 P) + 6 * penalty,
1475 lr[idx][j]);
1476 /* } */
1477 }
1478
1479 /* min_colonne=MIN2(lr[idx][j]+(type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]], min_colonne); */
1480 int bla;
1481 bla = lr[idx][j] + vrna_E_ext_stem(rtype[type], SS2[j - 1], SS1[i + 1], P) + 2 * penalty;
1482 min_colonne = MIN2(bla, min_colonne);
1483 }
1484 position[i] = min_colonne;
1485 if (max >= min_colonne) {
1486 max = min_colonne;
1487 max_pos = i;
1488 }
1489
1490 min_colonne = INF;
1491 }
1492
1493 free(S1);
1494 free(S2);
1495 free(SS1);
1496 free(SS2);
1497 if (max < threshTE) {
1498 find_max_snoop(s1,
1499 s2,
1500 max,
1501 alignment_length,
1502 position,
1503 delta,
1504 distance,
1505 penalty,
1506 threshloop,
1507 threshLE,
1508 threshRE,
1509 threshDE,
1510 threshTE,
1511 threshSE,
1512 threshD,
1513 half_stem,
1514 max_half_stem,
1515 min_s2,
1516 max_s2,
1517 min_s1,
1518 max_s1,
1519 min_d1,
1520 min_d2,
1521 name,
1522 fullStemEnergy);
1523 }
1524
1525 for (i = 1; i < 6; i++) {
1526 free(lc[i]);
1527 free(lr[i]);
1528 free(lpair[i]);
1529 }
1530 free(lc[0]);
1531 free(lr[0]);
1532 free(lpair[0]);
1533 free(lc);
1534 free(lr);
1535 free(lpair);
1536 free(position);
1537 }
1538
1539
1540 PRIVATE void
find_max_snoop(const char * s1,const char * s2,const int max,const int alignment_length,const int * position,const int delta,const int distance,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const char * name,const int fullStemEnergy)1541 find_max_snoop(const char *s1,
1542 const char *s2,
1543 const int max,
1544 const int alignment_length,
1545 const int *position,
1546 const int delta,
1547 const int distance,
1548 const int penalty,
1549 const int threshloop,
1550 const int threshLE,
1551 const int threshRE,
1552 const int threshDE,
1553 const int threshTE,
1554 const int threshSE,
1555 const int threshD,
1556 const int half_stem,
1557 const int max_half_stem,
1558 const int min_s2,
1559 const int max_s2,
1560 const int min_s1,
1561 const int max_s1,
1562 const int min_d1,
1563 const int min_d2,
1564 const char *name,
1565 const int fullStemEnergy)
1566 {
1567 int count = 0;
1568 int pos = n1 + 1;
1569 int threshold = MIN2(threshTE, max + delta);
1570
1571 /* printf("threshTE %d max %d\n", threshTE, max); */
1572 /* #pragma omp parallel for */
1573 /* for(pos=n1+1;pos>distance;pos--){ */
1574 while (pos-- > 5) {
1575 int temp_min = 0;
1576 if (position[pos] < (threshold)) {
1577 int search_range;
1578 search_range = distance + 1;
1579 while (--search_range)
1580 if (position[pos - search_range] <= position[pos - temp_min])
1581 temp_min = search_range;
1582
1583 pos -= temp_min;
1584 int begin = MAX2(6, pos - alignment_length + 1);
1585 char *s3 = (char *)vrna_alloc(sizeof(char) * (pos - begin + 3 + 12));
1586 strcpy(s3, "NNNNN");
1587 strncat(s3, (s1 + begin - 1), pos - begin + 2);
1588 strcat(s3, "NNNNN\0");
1589 /* printf("%s s3\n", s3); */
1590 snoopT test;
1591 test = snoopfold(s3,
1592 s2,
1593 penalty,
1594 threshloop,
1595 threshLE,
1596 threshRE,
1597 threshDE,
1598 threshD,
1599 half_stem,
1600 max_half_stem,
1601 min_s2,
1602 max_s2,
1603 min_s1,
1604 max_s1,
1605 min_d1,
1606 min_d2,
1607 fullStemEnergy);
1608 if (test.energy == INF) {
1609 free(s3);
1610 continue;
1611 }
1612
1613 if (test.Duplex_El > threshLE * 0.01 || test.Duplex_Er > threshRE * 0.01 ||
1614 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
1615 (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE * 0.01) {
1616 free(test.structure);
1617 free(s3);
1618 continue;
1619 }
1620
1621 int l1;
1622 l1 = strchr(test.structure, '&') - test.structure;
1623
1624 int shift = 0;
1625 if (test.i > (int)strlen(s3) - 10) {
1626 test.i--;
1627 l1--;
1628 }
1629
1630 if (test.i - l1 < 0) {
1631 l1--;
1632 shift++;
1633 }
1634
1635 char *target_struct = (char *)vrna_alloc(sizeof(char) * (strlen(test.structure) + 1));
1636 strncpy(target_struct, test.structure + shift, l1);
1637 strncat(target_struct, test.structure + (strchr(test.structure, '&') -
1638 test.structure),
1639 (int)strlen(test.structure) - (strchr(test.structure, '&') -
1640 test.
1641 structure));
1642 strcat(target_struct, "\0");
1643 char *target;
1644 target = (char *)vrna_alloc(l1 + 1);
1645 strncpy(target, (s3 + test.i + 5 - l1), l1);
1646 target[l1] = '\0';
1647 char *s4;
1648 s4 = (char *)vrna_alloc(sizeof(char) * (strlen(s2) - 9));
1649 strncpy(s4, s2 + 5, (int)strlen(s2) - 10);
1650 s4[(int)strlen(s2) - 10] = '\0';
1651 printf(
1652 "%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + 4.1 ) (%5.2f) \n%s&%s\n",
1653 target_struct,
1654 begin + test.i - 5 - l1,
1655 begin + test.i - 6,
1656 begin + test.u - 6,
1657 test.j + 1,
1658 test.j + (int)(strrchr(test.structure, '>') - strchr(test.structure, '>')) + 1,
1659 test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10,
1660 test.Duplex_El,
1661 test.Duplex_Er,
1662 test.Loop_E,
1663 test.Loop_D,
1664 test.fullStemEnergy,
1665 target,
1666 s4);
1667 if (name) {
1668 char *temp_seq;
1669 char *temp_struc;
1670 char *psoutput;
1671 temp_seq = (char *)vrna_alloc(sizeof(char) * (l1 + n2 - 9));
1672 temp_struc = (char *)vrna_alloc(sizeof(char) * (l1 + n2 - 9));
1673 strcpy(temp_seq, target);
1674 strcat(temp_seq, s4);
1675 strncpy(temp_struc, target_struct, l1);
1676 strcat(temp_struc, target_struct + l1 + 1);
1677 temp_seq[n2 + l1 - 10] = '\0';
1678 temp_struc[n2 + l1 - 10] = '\0';
1679 cut_point = l1 + 1;
1680 psoutput = vrna_strdup_printf("sno_%d_u_%d_%s.ps",
1681 count,
1682 begin + test.u - 6,
1683 name);
1684
1685 PS_rna_plot_snoop_a(temp_seq, temp_struc, psoutput, NULL, NULL);
1686 cut_point = -1;
1687 free(temp_seq);
1688 free(temp_struc);
1689 free(psoutput);
1690 count++;
1691 /* free(psoutput); */
1692 }
1693
1694 free(s4);
1695 free(test.structure);
1696 free(target_struct);
1697 free(target);
1698 free(s3);
1699 }
1700 }
1701 }
1702
1703
1704 snoopT
snoopfold(const char * s1,const char * s2,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const int fullStemEnergy)1705 snoopfold(const char *s1,
1706 const char *s2,
1707 const int penalty,
1708 const int threshloop,
1709 const int threshLE,
1710 const int threshRE,
1711 const int threshDE,
1712 const int threshD,
1713 const int half_stem,
1714 const int max_half_stem,
1715 const int min_s2,
1716 const int max_s2,
1717 const int min_s1,
1718 const int max_s1,
1719 const int min_d1,
1720 const int min_d2,
1721 const int fullStemEnergy)
1722 {
1723 /* int Eminj, Emin_l; */
1724 int i, j, l1, Emin = INF, i_min = 0, j_min = 0;
1725 char *struc;
1726 snoopT mfe;
1727 int *indx;
1728 int *mLoop;
1729 int *cLoop;
1730 folden **foldlist, **foldlist_XS;
1731 int Duplex_El, Duplex_Er;
1732 int Loop_D;
1733 int u;
1734 int Loop_E;
1735 vrna_md_t md;
1736
1737 Duplex_El = 0;
1738 Duplex_Er = 0;
1739 Loop_E = 0, Loop_D = 0;
1740 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
1741 n1 = (int)strlen(s1);
1742 n2 = (int)strlen(s2);
1743
1744 set_model_details(&md);
1745 if ((!P) || (fabs(P->temperature - temperature) > 1e-6)) {
1746 snoupdate_fold_params();
1747 if (P)
1748 free(P);
1749
1750 P = vrna_params(&md);
1751 make_pair_matrix();
1752 }
1753
1754 c = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
1755 r = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
1756 for (i = 0; i <= n1; i++) {
1757 c[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1758 r[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1759 for (j = n2; j > -1; j--) {
1760 c[i][j] = INF;
1761 r[i][j] = INF;
1762 }
1763 }
1764 encode_seqs(s1, s2);
1765 for (i = 6; i <= n1 - 5; i++) {
1766 for (j = n2 - min_d2; j > min_d1; j--) {
1767 int type, type2, E, k, l;
1768 type = pair[S1[i]][S2[j]];
1769 c[i][j] = (type) ? P->DuplexInit : INF;
1770 if (!type)
1771 continue;
1772
1773 if (/* pair[S1[i+1]][S2[j-1]] && */
1774 j < max_s1 && j > min_s1 &&
1775 j > n2 - max_s2 - max_half_stem &&
1776 j < n2 - min_s2 - half_stem && S1[i - 2] == 4) {
1777 /*constraint on s2 and i*/
1778 int min_k, max_k;
1779 max_k = MIN2(n2 - min_s2, j + max_half_stem);
1780 min_k = MAX2(j + half_stem, n2 - max_s2);
1781 folden *temp;
1782 temp = foldlist[j + 1];
1783 while (temp->next) {
1784 int k = temp->k;
1785 /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
1786 if (pair[S1[i - 3]][S2[k + 1]] /*&& pair[S1[i-4]][S2[k+2]]*/)
1787 r[i][j] = MIN2(r[i][j], c[i - 3][k + 1] + temp->energy);
1788
1789 /*else*/ if (pair[S1[i - 4]][S2[k + 1]] /*&& pair[S1[i-5]][S2[k+2]]*/)
1790 r[i][j] = MIN2(r[i][j], c[i - 4][k + 1] + temp->energy);
1791
1792 /* } */
1793 temp = temp->next;
1794 }
1795 }
1796
1797 /* dangle 5'SIDE relative to the mRNA */
1798 /**
1799 *** c[i][j] += P->dangle5[type][SS1[i-1]];
1800 *** c[i][j] += P->dangle3[type][SS2[j+1]];
1801 *** if (type>2) c[i][j] += P->TerminalAU;
1802 **/
1803 c[i][j] += vrna_E_ext_stem(type, SS1[i - 1], SS2[j + 1], P);
1804 for (k = i - 1; k > 0 && (i - k) < MAXLOOP_L; k--) {
1805 for (l = j + 1; l <= n2; l++) {
1806 if (i - k + l - j > 2 * MAXLOOP_L - 2)
1807 break;
1808
1809 if (abs(i - k - l + j) >= ASS)
1810 continue;
1811
1812 type2 = pair[S1[k]][S2[l]];
1813 if (!type2)
1814 continue;
1815
1816 E = E_IntLoop(i - k - 1, l - j - 1, type2, rtype[type],
1817 SS1[k + 1], SS2[l - 1], SS1[i - 1], SS2[j + 1], P);
1818 c[i][j] = MIN2(c[i][j], c[k][l] + E + (i - k + l - j) * penalty);
1819 r[i][j] = MIN2(r[i][j], r[k][l] + E + (i - k + l - j) * penalty);
1820 }
1821 }
1822 E = r[i][j];
1823 /**
1824 *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
1825 *** if (j>1) E += P->dangle5[rtype[type]][SS2[j-1]];
1826 *** f (type>2) E += P->TerminalAU;
1827 **/
1828 E += vrna_E_ext_stem(rtype[type], (j > 1) ? SS2[j - 1] : -1, (i < n1) ? SS1[i + 1] : -1, P);
1829 if (E < Emin) {
1830 Emin = E;
1831 i_min = i;
1832 j_min = j;
1833 }
1834 }
1835 }
1836 if (Emin > 0) {
1837 printf("no target found under the constraints chosen\n");
1838 for (i = 0; i <= n1; i++) {
1839 free(r[i]);
1840 free(c[i]);
1841 }
1842 free(c);
1843 free(r);
1844 free(S1);
1845 free(S2);
1846 free(SS1);
1847 free(SS2);
1848 mfe.energy = INF;
1849 return mfe;
1850 }
1851
1852 struc = snoop_backtrack(i_min, j_min, s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D,
1853 &u, penalty, threshloop, threshLE, threshRE, threshDE, threshD,
1854 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
1855 /* if (i_min<n1-5) i_min++; */
1856 /* if (j_min>1 ) j_min--; */
1857 l1 = strchr(struc, '&') - struc;
1858 mfe.i = i_min - 5;
1859 mfe.j = j_min - 5;
1860 mfe.u = u - 5;
1861 mfe.Duplex_Er = (float)Duplex_Er / 100;
1862 mfe.Duplex_El = (float)Duplex_El / 100;
1863 mfe.Loop_D = (float)Loop_D / 100;
1864 mfe.Loop_E = (float)Loop_E / 100;
1865 mfe.energy = (float)Emin / 100;
1866 mfe.fullStemEnergy = (float)fullStemEnergy / 100;
1867 mfe.structure = struc;
1868 if (!delay_free) {
1869 for (i = 0; i <= n1; i++) {
1870 free(r[i]);
1871 free(c[i]);
1872 }
1873 free(c);
1874 free(r);
1875 free(S1);
1876 free(S2);
1877 free(SS1);
1878 free(SS2);
1879 }
1880
1881 return mfe;
1882 }
1883
1884
1885 PRIVATE int
snoopfold_XS_fill(const char * s1,const char * s2,const int ** access_s1,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2)1886 snoopfold_XS_fill(const char *s1,
1887 const char *s2,
1888 const int **access_s1,
1889 const int penalty,
1890 const int threshloop,
1891 const int threshLE,
1892 const int threshRE,
1893 const int threshDE,
1894 const int threshD,
1895 const int half_stem,
1896 const int max_half_stem,
1897 const int min_s2,
1898 const int max_s2,
1899 const int min_s1,
1900 const int max_s1,
1901 const int min_d1,
1902 const int min_d2)
1903 {
1904 /* int Eminj, Emin_l; */
1905 int i, j, Emin = INF, i_min = 0, j_min = 0;
1906 /* char *struc; */
1907 /* snoopT mfe; */
1908 int *indx;
1909 int *mLoop;
1910 int *cLoop;
1911 folden **foldlist, **foldlist_XS;
1912 int Duplex_El, Duplex_Er;
1913 int Loop_D;
1914 /* int u; */
1915 int Loop_E;
1916 vrna_md_t md;
1917
1918 Duplex_El = 0;
1919 Duplex_Er = 0;
1920 Loop_E = 0, Loop_D = 0;
1921 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
1922 n1 = (int)strlen(s1);
1923 n2 = (int)strlen(s2);
1924
1925 set_model_details(&md);
1926 if ((!P) || (fabs(P->temperature - temperature) > 1e-6)) {
1927 snoupdate_fold_params();
1928 if (P)
1929 free(P);
1930
1931 P = vrna_params(&md);
1932 make_pair_matrix();
1933 }
1934
1935 c_fill = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
1936 r_fill = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
1937 for (i = 0; i <= n1; i++) {
1938 c_fill[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1939 r_fill[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
1940 for (j = n2; j > -1; j--) {
1941 c_fill[i][j] = INF;
1942 r_fill[i][j] = INF;
1943 }
1944 }
1945 encode_seqs(s1, s2);
1946
1947 int di[5];
1948 di[0] = 0;
1949 for (i = 6; i <= n1 - 5; i++) {
1950 di[1] = access_s1[5][i] - access_s1[4][i - 1];
1951 di[2] = access_s1[5][i - 1] - access_s1[4][i - 2] + di[1];
1952 di[3] = access_s1[5][i - 2] - access_s1[4][i - 3] + di[2];
1953 di[4] = access_s1[5][i - 3] - access_s1[4][i - 4] + di[3];
1954 di[1] = MIN2(di[1], 165);
1955 di[2] = MIN2(di[2], 330);
1956 di[3] = MIN2(di[3], 495);
1957 di[4] = MIN2(di[4], 660);
1958 for (j = n2 - min_d2; j > min_d1; j--) {
1959 int type, type2, E, k, l;
1960 type = pair[S1[i]][S2[j]];
1961 c_fill[i][j] = (type) ? P->DuplexInit : INF;
1962 if (!type)
1963 continue;
1964
1965 if (/* pair[S1[i+1]][S2[j-1]] && */
1966 j < max_s1 && j > min_s1 &&
1967 j > n2 - max_s2 - max_half_stem &&
1968 j < n2 - min_s2 - half_stem && S1[i - 2] == 4) {
1969 /*constraint on s2 and i*/
1970 int min_k, max_k;
1971 max_k = MIN2(n2 - min_s2, j + max_half_stem);
1972 min_k = MAX2(j + half_stem, n2 - max_s2);
1973 folden *temp;
1974 temp = foldlist[j + 1];
1975 while (temp->next) {
1976 int k = temp->k;
1977 /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
1978 if (pair[S1[i - 3]][S2[k + 1]] /*&& pair[S1[i-4]][S2[k+2]]*/)
1979 r_fill[i][j] = MIN2(r_fill[i][j], c_fill[i - 3][k + 1] + temp->energy + di[3]);
1980
1981 /*else*/ if (pair[S1[i - 4]][S2[k + 1]] /*&& pair[S1[i-5]][S2[k+2]]*/)
1982 r_fill[i][j] = MIN2(r_fill[i][j], c_fill[i - 4][k + 1] + temp->energy + di[4]);
1983
1984 /* } */
1985 temp = temp->next;
1986 }
1987 }
1988
1989 /* dangle 5'SIDE relative to the mRNA */
1990 /**
1991 *** c_fill[i][j] += P->dangle5[type][SS1[i-1]];
1992 *** c_fill[i][j] += P->dangle3[type][SS2[j+1]];
1993 *** if (type>2) c_fill[i][j] += P->TerminalAU;
1994 **/
1995 c_fill[i][j] += vrna_E_ext_stem(type, SS1[i - 1], SS2[j + 1], P);
1996 for (k = i - 1; k > 0 && (i - k) < MAXLOOP_L; k--) {
1997 for (l = j + 1; l <= n2; l++) {
1998 if (i - k + l - j > 2 * MAXLOOP_L - 2)
1999 break;
2000
2001 if (abs(i - k - l + j) >= ASS)
2002 continue;
2003
2004 type2 = pair[S1[k]][S2[l]];
2005 if (!type2)
2006 continue;
2007
2008 E = E_IntLoop(i - k - 1, l - j - 1, type2, rtype[type],
2009 SS1[k + 1], SS2[l - 1], SS1[i - 1], SS2[j + 1], P);
2010 c_fill[i][j] = MIN2(c_fill[i][j], c_fill[k][l] + E + di[i - k]);
2011 r_fill[i][j] = MIN2(r_fill[i][j], r_fill[k][l] + E + di[i - k]);
2012 }
2013 }
2014 E = r_fill[i][j];
2015 /**
2016 *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
2017 *** if (j>1) E += P->dangle5[rtype[type]][SS2[j-1]];
2018 *** if (type>2) E += P->TerminalAU;
2019 **/
2020 E += vrna_E_ext_stem(rtype[type], (j > 1) ? SS2[j - 1] : -1, (i < n1) ? SS1[i + 1] : -1, P);
2021 if (E < Emin) {
2022 Emin = E;
2023 i_min = i;
2024 j_min = j;
2025 }
2026 }
2027 }
2028 return Emin;
2029 }
2030
2031
2032 PUBLIC snoopT *
snoop_subopt(const char * s1,const char * s2,int delta,int w,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int distance,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const int fullStemEnergy)2033 snoop_subopt(const char *s1,
2034 const char *s2,
2035 int delta,
2036 int w,
2037 const int penalty,
2038 const int threshloop,
2039 const int threshLE,
2040 const int threshRE,
2041 const int threshDE,
2042 const int threshTE,
2043 const int threshSE,
2044 const int threshD,
2045 const int distance,
2046 const int half_stem,
2047 const int max_half_stem,
2048 const int min_s2,
2049 const int max_s2,
2050 const int min_s1,
2051 const int max_s1,
2052 const int min_d1,
2053 const int min_d2,
2054 const int fullStemEnergy)
2055 {
2056 /* printf("%d %d\n", min_s2, max_s2); */
2057 int i, j, n1, n2, E, n_subopt = 0, n_max;
2058 char *struc;
2059 snoopT mfe;
2060 snoopT *subopt;
2061 int thresh;
2062
2063 int Duplex_El, Duplex_Er, Loop_E;
2064 int Loop_D;
2065
2066 Duplex_El = 0;
2067 Duplex_Er = 0;
2068 Loop_E = 0;
2069 Loop_D = 0;
2070 int u;
2071 u = 0;
2072 n_max = 16;
2073 subopt = (snoopT *)vrna_alloc(n_max * sizeof(snoopT));
2074 delay_free = 1;
2075 mfe = snoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE, threshD,
2076 half_stem, max_half_stem,
2077 min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, fullStemEnergy);
2078
2079
2080 if (mfe.energy > 0) {
2081 free(subopt);
2082 delay_free = 0;
2083 return NULL;
2084 }
2085
2086 thresh = MIN2((int)((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E) * 100 + 0.1 + 410) + delta,
2087 threshTE);
2088 /* subopt[n_subopt++]=mfe; */
2089 free(mfe.structure);
2090
2091 n1 = (int)strlen(s1);
2092 n2 = (int)strlen(s2);
2093 for (i = n1; i > 0; i--) {
2094 for (j = 1; j <= n2; j++) {
2095 int type, Ed;
2096 type = pair[S2[j]][S1[i]];
2097 if (!type)
2098 continue;
2099
2100 E = Ed = r[i][j];
2101 /**
2102 *** if (i<n1) Ed += P->dangle3[type][SS1[i+1]];
2103 *** if (j>1) Ed += P->dangle5[type][SS2[j-1]];
2104 *** if (type>2) Ed += P->TerminalAU;
2105 **/
2106 Ed += vrna_E_ext_stem(type, (j > 1) ? SS2[j - 1] : -1, (i < n1) ? SS1[i + 1] : -1, P);
2107 if (Ed > thresh)
2108 continue;
2109
2110 /* too keep output small, remove hits that are dominated by a
2111 * better one close (w) by. For simplicity we do test without
2112 * adding dangles, which is slightly inaccurate.
2113 */
2114 /* w=1; */
2115 /* for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) { */
2116 /* for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++) */
2117 /* if (r[ii][jj]<E) {type=0; break;} */
2118 /* } */
2119 if (!type)
2120 continue;
2121
2122 struc = snoop_backtrack(i,
2123 j,
2124 s2,
2125 &Duplex_El,
2126 &Duplex_Er,
2127 &Loop_E,
2128 &Loop_D,
2129 &u,
2130 penalty,
2131 threshloop,
2132 threshLE,
2133 threshRE,
2134 threshDE,
2135 threshD,
2136 half_stem,
2137 max_half_stem,
2138 min_s2,
2139 max_s2,
2140 min_s1,
2141 max_s1,
2142 min_d1,
2143 min_d2);
2144 if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
2145 (Duplex_Er + Duplex_El) > threshDE ||
2146 (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
2147 (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
2148 /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
2149 /* " Duplex_Er + Duplex_El %d threshDE %d \n" */
2150 /* " Duplex_Er + Duplex_El + Loop_E %d threshTE %d \n" */
2151 /* " Duplex_Er + Duplex_El + Loop_E + Loop_D %d threshSE %d \n", */
2152 /* Duplex_Er , threshRE , Duplex_El ,threshLE, */
2153 /* Duplex_Er + Duplex_El, threshDE, */
2154 /* Duplex_Er + Duplex_El+ Loop_E , threshTE, */
2155 /* Duplex_Er + Duplex_El+ Loop_E + Loop_D, threshSE); */
2156 Duplex_Er = 0;
2157 Duplex_El = 0;
2158 Loop_E = 0;
2159 Loop_D = 0;
2160 u = 0,
2161 free(struc);
2162 continue;
2163 }
2164
2165 if (n_subopt + 1 >= n_max) {
2166 n_max *= 2;
2167 subopt = (snoopT *)vrna_realloc(subopt, n_max * sizeof(snoopT));
2168 }
2169
2170 subopt[n_subopt].i = i - 5;
2171 subopt[n_subopt].j = j - 5;
2172 subopt[n_subopt].u = u - 5;
2173 subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
2174 subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
2175 subopt[n_subopt].Loop_E = Loop_E * 0.01;
2176 subopt[n_subopt].Loop_D = Loop_D * 0.01;
2177 subopt[n_subopt].energy = (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) * 0.01;
2178 subopt[n_subopt].fullStemEnergy = (float)fullStemEnergy * 0.01;
2179 subopt[n_subopt++].structure = struc;
2180
2181 Duplex_Er = 0;
2182 Duplex_El = 0;
2183 Loop_E = 0;
2184 Loop_D = 0;
2185 u = 0;
2186 }
2187 }
2188
2189 for (i = 0; i <= n1; i++) {
2190 free(c[i]);
2191 free(r[i]);
2192 }
2193 free(c);
2194 free(r);
2195 free(S1);
2196 free(S2);
2197 free(SS1);
2198 free(SS2);
2199 delay_free = 0;
2200
2201 if (snoop_subopt_sorted)
2202 qsort(subopt, n_subopt, sizeof(snoopT), compare);
2203
2204 subopt[n_subopt].i = 0;
2205 subopt[n_subopt].j = 0;
2206 subopt[n_subopt].structure = NULL;
2207 return subopt;
2208 }
2209
2210
2211 PUBLIC void
snoop_subopt_XS(const char * s1,const char * s2,const int ** access_s1,int delta,int w,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int distance,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const int alignment_length,const char * name,const int fullStemEnergy)2212 snoop_subopt_XS(const char *s1,
2213 const char *s2,
2214 const int **access_s1,
2215 int delta,
2216 int w,
2217 const int penalty,
2218 const int threshloop,
2219 const int threshLE,
2220 const int threshRE,
2221 const int threshDE,
2222 const int threshTE,
2223 const int threshSE,
2224 const int threshD,
2225 const int distance,
2226 const int half_stem,
2227 const int max_half_stem,
2228 const int min_s2,
2229 const int max_s2,
2230 const int min_s1,
2231 const int max_s1,
2232 const int min_d1,
2233 const int min_d2,
2234 const int alignment_length,
2235 const char *name,
2236 const int fullStemEnergy)
2237 {
2238 /* printf("%d %d\n", min_s2, max_s2); */
2239 int i, j, E, n_max;
2240 /* char *struc; */
2241 /* snoopT mfe; */
2242
2243 int thresh;
2244
2245 int Duplex_El, Duplex_Er, Loop_E;
2246 int Loop_D;
2247
2248 Duplex_El = 0;
2249 Duplex_Er = 0;
2250 Loop_E = 0;
2251 Loop_D = 0;
2252 int u;
2253 u = 0;
2254 n_max = 16;
2255 delay_free = 1;
2256 int Emin = snoopfold_XS_fill(s1,
2257 s2,
2258 access_s1,
2259 penalty,
2260 threshloop,
2261 threshLE,
2262 threshRE,
2263 threshDE,
2264 threshD,
2265 half_stem,
2266 max_half_stem,
2267 min_s2,
2268 max_s2,
2269 min_s1,
2270 max_s1,
2271 min_d1,
2272 min_d2);
2273 if (Emin > 0)
2274 delay_free = 0;
2275
2276 thresh = MIN2(-100, threshTE + alignment_length * 30);
2277 /* n1=(int)strlen(s1); */
2278 /* n2=(int)strlen(s2); */
2279
2280 int n3 = (int)strlen(s1);
2281 int n4 = (int)strlen(s2);
2282 S1_fill = (short *)vrna_alloc(sizeof(short) * (n3 + 2));
2283 S2_fill = (short *)vrna_alloc(sizeof(short) * (n4 + 2));
2284 SS1_fill = (short *)vrna_alloc(sizeof(short) * (n3 + 1));
2285 SS2_fill = (short *)vrna_alloc(sizeof(short) * (n4 + 1));
2286 memcpy(S1_fill, S1, sizeof(short) * n3 + 2);
2287 memcpy(S2_fill, S2, sizeof(short) * n4 + 2);
2288 memcpy(SS1_fill, SS1, sizeof(short) * n3 + 1);
2289 memcpy(SS2_fill, SS2, sizeof(short) * n4 + 1);
2290 free(S1);
2291 free(S2);
2292 free(SS1);
2293 free(SS2);
2294 int count = 0;
2295 for (i = n3 - 5; i > 0; i--) {
2296 for (j = 1; j <= n4; j++) {
2297 int type, Ed;
2298 type = pair[S2_fill[j]][S1_fill[i]];
2299 if (!type)
2300 continue;
2301
2302 E = Ed = r_fill[i][j];
2303 /**
2304 ***if (i<n3) Ed += P->dangle3[type][SS1_fill[i+1]];
2305 ***if (j>1) Ed += P->dangle5[type][SS2_fill[j-1]];
2306 ***if (type>2) Ed += P->TerminalAU;
2307 **/
2308 Ed += vrna_E_ext_stem(type, (j > 1) ? SS2[j - 1] : -1, (i < n3) ? SS1[i + 1] : -1, P);
2309 if (Ed > thresh)
2310 continue;
2311
2312 /* to keep output small, remove hits that are dominated by a
2313 * better one close (w) by. For simplicity we do test without
2314 * adding dangles, which is slightly inaccurate.
2315 */
2316 /* w=10; */
2317 /* for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n3-5)) && type; ii++) { */
2318 /* for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n4-5); jj++) */
2319 /* if (r_fill[ii][jj]<E) {type=0; break;} */
2320 /* } */
2321 /* i=ii;j=jj; */
2322 if (!type)
2323 continue;
2324
2325 int begin = MAX2(5, i - alignment_length);
2326 int end = MIN2(n3 - 5, i - 1);
2327 char *s3 = (char *)vrna_alloc(sizeof(char) * (end - begin + 2) + 5);
2328 strncpy(s3, (s1 + begin), end - begin + 1);
2329 strcat(s3, "NNNNN\0");
2330 int n5 = (int)strlen(s3);
2331 snoopT test = snoopfold_XS(s3, s2, access_s1, i, j, penalty,
2332 threshloop, threshLE, threshRE,
2333 threshDE, threshD, half_stem,
2334 max_half_stem, min_s2, max_s2, min_s1,
2335 max_s1, min_d1, min_d2, fullStemEnergy);
2336 if (test.energy == INF) {
2337 free(s3);
2338 continue;
2339 }
2340
2341 if (test.Duplex_El > threshLE * 0.01 || test.Duplex_Er > threshRE * 0.01 ||
2342 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
2343 (test.Duplex_Er + test.Duplex_El + test.Loop_E) > threshTE * 0.01 ||
2344 (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE * 0.01) {
2345 free(test.structure);
2346 free(s3);
2347 continue;
2348 }
2349
2350 char *s4;
2351 s4 = (char *)vrna_alloc(sizeof(char) * (n4 - 9));
2352 strncpy(s4, s2 + 5, n4 - 10);
2353 s4[n4 - 10] = '\0';
2354
2355 char *s5 = vrna_alloc(sizeof(char) * n5 - test.i + 2 - 5);
2356 strncpy(s5, s3 + test.i - 1, n5 - test.i + 1 - 5);
2357 s5[n5 - test.i + 1 - 5] = '\0';
2358 float dE = ((float)(access_s1[n5 - test.i + 1 - 5][i])) * 0.01;
2359 printf(
2360 "%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10) (%5.2f)\n%s&%s\n",
2361 test.structure,
2362 i - (n5 - test.i),
2363 i - 5,
2364 i - (n5 - test.u),
2365 j - 5,
2366 j - 5 + (int)(strrchr(test.structure, '>') - strchr(test.structure, '>')),
2367 test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10 + dE,
2368 test.Duplex_El,
2369 test.Duplex_Er,
2370 test.Loop_E,
2371 test.Loop_D,
2372 dE,
2373 test.fullStemEnergy,
2374 s5,
2375 s4);
2376 if (name) {
2377 int begin_t, end_t, begin_q, end_q, and, pipe, k;
2378 char *psoutput;
2379 begin_q = 0;
2380 end_q = n4 - 10;
2381 begin_t = 0;
2382 end_t = n5 - test.i + 1 - 5;
2383 and = end_t + 1;
2384 pipe = test.u - test.i + 1;
2385 cut_point = end_t + 1;
2386 char *catseq, *catstruct;/* *fname; */
2387 catseq = (char *)vrna_alloc(n5 + end_q - begin_q + 2);
2388 catstruct = (char *)vrna_alloc(n5 + end_q - begin_q + 2);
2389 strcpy(catseq, s5);
2390 strncpy(catstruct, test.structure, end_t);
2391 strcat(catseq, s4);
2392 strncat(catstruct, test.structure + end_t + 1, end_q - begin_q + 1);
2393 catstruct[end_t - begin_t + end_q - begin_q + 2] = '\0';
2394 catseq[end_t - begin_t + end_q - begin_q + 2] = '\0';
2395 int *relative_access;
2396 relative_access = vrna_alloc(sizeof(int) * strlen(s5));
2397 relative_access[0] = access_s1[1][i - (n5 - test.i) + 5];
2398 for (k = 1; k < (int)strlen(s5); k++)
2399 relative_access[k] = access_s1[k + 1][i - (n5 - test.i) + k + 5] -
2400 access_s1[k][i - (n5 - test.i) + k + 4];
2401
2402 psoutput = vrna_strdup_printf("sno_XS_%d_u_%d_%s.ps",
2403 count,
2404 i - (n5 - test.u),
2405 name);
2406
2407 PS_rna_plot_snoop_a(catseq, catstruct, psoutput, relative_access, NULL);
2408 free(catseq);
2409 free(catstruct);
2410 free(relative_access);
2411 free(psoutput);
2412 count++;
2413 }
2414
2415 free(s3);
2416 free(s4);
2417 free(s5);
2418 free(test.structure);
2419 }
2420 }
2421 for (i = 0; i <= n3; i++) {
2422 free(c_fill[i]);
2423 free(r_fill[i]);
2424 }
2425 free(c_fill);
2426 free(r_fill);
2427 free(S1_fill);
2428 free(S2_fill);
2429 free(SS1_fill);
2430 free(SS2_fill);
2431 delay_free = 0;
2432 }
2433
2434
2435 PRIVATE char *
snoop_backtrack(int i,int j,const char * snoseq,int * Duplex_El,int * Duplex_Er,int * Loop_E,int * Loop_D,int * u,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2)2436 snoop_backtrack(int i,
2437 int j,
2438 const char *snoseq,
2439 int *Duplex_El,
2440 int *Duplex_Er,
2441 int *Loop_E,
2442 int *Loop_D,
2443 int *u,
2444 const int penalty,
2445 const int threshloop,
2446 const int threshLE,
2447 const int threshRE,
2448 const int threshDE,
2449 const int threshD,
2450 const int half_stem,
2451 const int max_half_stem,
2452 const int min_s2,
2453 const int max_s2,
2454 const int min_s1,
2455 const int max_s1,
2456 const int min_d1,
2457 const int min_d2)
2458 {
2459 /* backtrack structure going backwards from i, and forwards from j
2460 * return structure in bracket notation with & as separator */
2461 int k, l, type, type2, E, traced, i0, j0;
2462 int traced_r = 0; /* flag for following backtrack in c or r */
2463 char *st1, *st2, *struc;
2464 char *struc_loop;
2465
2466 st1 = (char *)vrna_alloc(sizeof(char) * (n1 + 1));
2467 st2 = (char *)vrna_alloc(sizeof(char) * (n2 + 1));
2468 int *indx;
2469 int *mLoop;
2470 int *cLoop;
2471 folden **foldlist, **foldlist_XS;
2472 type = pair[S1[i]][S2[j]];
2473 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
2474 i0 = i;
2475 j0 = j;
2476 /**
2477 *** if (i<n1) *Duplex_Er += P->dangle3[rtype[type]][SS1[i+1]];
2478 *** if (j>1) *Duplex_Er += P->dangle5[rtype[type]][SS2[j-1]];
2479 *** if (type>2) *Duplex_Er += P->TerminalAU;
2480 **/
2481 *Duplex_Er += vrna_E_ext_stem(rtype[type], (j > 1) ? SS2[j - 1] : -1, (i < n1) ? SS1[i + 1] : -1, P);
2482 while (i > 0 && j <= n2 - min_d2) {
2483 if (!traced_r) {
2484 E = r[i][j];
2485 traced = 0;
2486 st1[i - 1] = '<';
2487 st2[j - 1] = '>';
2488 type = pair[S1[i]][S2[j]];
2489 if (!type)
2490 vrna_message_error("backtrack failed in fold duplex r");
2491
2492 for (k = i - 1; k > 0 && (i - k) < MAXLOOP_L; k--) {
2493 for (l = j + 1; l <= n2; l++) {
2494 int LE;
2495 if (i - k + l - j > 2 * MAXLOOP_L - 2)
2496 break;
2497
2498 if (abs(i - k - l + j) >= ASS)
2499 continue;
2500
2501 type2 = pair[S1[k]][S2[l]];
2502 if (!type2)
2503 continue;
2504
2505 LE = E_IntLoop(i - k - 1, l - j - 1, type2, rtype[type],
2506 SS1[k + 1], SS2[l - 1], SS1[i - 1], SS2[j + 1], P);
2507 if (E == r[k][l] + LE + (i - k + l - j) * penalty) {
2508 traced = 1;
2509 i = k;
2510 j = l;
2511 *Duplex_Er += LE;
2512 break;
2513 }
2514 }
2515 if (traced)
2516 break;
2517 }
2518 if (!traced) {
2519 if (/* pair[S1[i+1]][S2[j-1]] && */
2520 j < max_s1 && j > min_s1 &&
2521 j > n2 - max_s2 - max_half_stem &&
2522 j < n2 - min_s2 - half_stem &&
2523 S1[i - 2] == 4) {
2524 int min_k, max_k;
2525 max_k = MIN2(n2 - min_s2, j + max_half_stem + 1);
2526 min_k = MAX2(j + half_stem + 1, n2 - max_s2);
2527 folden *temp;
2528 temp = foldlist[j + 1];
2529 while (temp->next) {
2530 int k = temp->k;
2531 if (pair[S1[i - 3]][S2[k + 1]] /*&& pair[S1[i-4]][S2[k+2]]*/) {
2532 /* introduce structure from RNAfold */
2533 if (E == c[i - 3][k + 1] + temp->energy) {
2534 *Loop_E = temp->energy;
2535 st1[i - 3] = '|';
2536 *u = i - 2;
2537 int a, b;
2538 /* int fix_ij=indx[k-1+1]+j+1; */
2539 for (a = 0; a < MISMATCH; a++) {
2540 for (b = 0; b < MISMATCH; b++) {
2541 int ij = indx[k - 1 - a + 1] + j + 1 + b;
2542 if (cLoop[ij] == temp->energy) {
2543 struc_loop = snobacktrack_fold_from_pair(snoseq, j + 1 + b, k - a - 1 + 1);
2544 a = INF;
2545 b = INF;
2546 }
2547 }
2548 }
2549 traced = 1;
2550 traced_r = 1;
2551 i = i - 3;
2552 j = k + 1;
2553 break;
2554 }
2555 }
2556
2557 /*else*/ if (pair[S1[i - 4]][S2[k + 1]] /*&& pair[S1[i-5]][S2[k+2]]*/) {
2558 /* introduce structure from RNAfold */
2559 if (E == c[i - 4][k + 1] + temp->energy) {
2560 *Loop_E = temp->energy;
2561 st1[i - 3] = '|';
2562 *u = i - 2;
2563 int a, b;
2564 /* int fix_ij=indx[k-1+1]+j+1; */
2565 for (a = 0; a < MISMATCH; a++) {
2566 for (b = 0; b < MISMATCH; b++) {
2567 int ij = indx[k - 1 - a + 1] + j + 1 + b;
2568 if (cLoop[ij] == temp->energy) {
2569 struc_loop = snobacktrack_fold_from_pair(snoseq, j + 1 + b, k - a - 1 + 1);
2570 a = INF;
2571 b = INF;
2572 }
2573 }
2574 }
2575 traced = 1;
2576 traced_r = 1;
2577 i = i - 4;
2578 j = k + 1;
2579 break;
2580 }
2581 } /* else if */
2582
2583 temp = temp->next;
2584 } /* while temp-> next */
2585 } /* test on j */
2586 } /* traced? */
2587 } /* traced_r? */
2588 else {
2589 E = c[i][j];
2590 traced = 0;
2591 st1[i - 1] = '<';
2592 st2[j - 1] = '>';
2593 type = pair[S1[i]][S2[j]];
2594 if (!type)
2595 vrna_message_error("backtrack failed in fold duplex c");
2596
2597 for (k = i - 1; (i - k) < MAXLOOP_L; k--) {
2598 for (l = j + 1; l <= n2; l++) {
2599 int LE;
2600 if (i - k + l - j > 2 * MAXLOOP_L - 2)
2601 break;
2602
2603 if (abs(i - k - l + j) >= ASS)
2604 continue;
2605
2606 type2 = pair[S1[k]][S2[l]];
2607 if (!type2)
2608 continue;
2609
2610 LE = E_IntLoop(i - k - 1, l - j - 1, type2, rtype[type],
2611 SS1[k + 1], SS2[l - 1], SS1[i - 1], SS2[j + 1], P);
2612 if (E == c[k][l] + LE + (i - k + l - j) * penalty) {
2613 traced = 1;
2614 i = k;
2615 j = l;
2616 *Duplex_El += LE;
2617 break;
2618 }
2619 }
2620 if (traced)
2621 break;
2622 }
2623 }
2624
2625 if (!traced) {
2626 int correction;
2627 correction = vrna_E_ext_stem(type, (i > 1) ? SS1[i - 1] : -1, (j < n2) ? SS2[j + 1] : -1, P);
2628 E -= correction;
2629 *Duplex_El += correction;
2630 /**
2631 *** if (i>1) {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];}
2632 *** if (j<n2) {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];}
2633 *** if (type>2) {E -= P->TerminalAU; *Duplex_El +=P->TerminalAU;}
2634 **/
2635 if (E != P->DuplexInit)
2636 vrna_message_error("backtrack failed in fold duplex end");
2637 else
2638 break;
2639 }
2640 }
2641 /* if (i>1) i--; */
2642 /* if (j<n2) j++; */
2643 /* struc = (char *) vrna_alloc(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
2644 struc = (char *)vrna_alloc(i0 - i + 1 + n2 - 1 + 1 + 2); /* declare final duplex structure */
2645 char *struc2;
2646 struc2 = (char *)vrna_alloc(n2 + 1);
2647 /* char * struct_const; */
2648 for (k = MAX2(i, 1); k <= i0; k++)
2649 if (!st1[k - 1])
2650 st1[k - 1] = '.';
2651
2652 /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
2653 /* char * struct_const; */
2654 /* struct_const = (char *) vrna_alloc(sizeof(char)*(n2+1)); */
2655 for (k = 1; k <= n2; k++) {
2656 if (!st2[k - 1])
2657 st2[k - 1] = struc_loop[k - 1]; /* '.'; */
2658
2659 struc2[k - 1] = st2[k - 1]; /* '.'; */
2660 /* if (k>=j0 && k<=j){ */
2661 /* struct_const[k-1]='x'; */
2662 /* } */
2663 /* else{ */
2664 /* if(k<j0) {struct_const[k-1]='<';} */
2665 /* if(k>j) {struct_const[k-1]='>';} */
2666 /* } */
2667 }
2668 char duplexseq_1[j0];
2669 char duplexseq_2[n2 - j + 2];
2670 if (j < n2) {
2671 strncpy(duplexseq_1, snoseq, j0 - 1);
2672 strcpy(duplexseq_2, snoseq + j);
2673 duplexseq_1[j0 - 1] = '\0';
2674 duplexseq_2[n2 - j + 1] = '\0';
2675 duplexT temp;
2676 temp = duplexfold(duplexseq_1, duplexseq_2);
2677 *Loop_D = MIN2(0, -410 + (int)100 * temp.energy);
2678 if (*Loop_D) {
2679 int l1, ibegin, iend, jbegin, jend;
2680 l1 = strchr(temp.structure, '&') - temp.structure;
2681 ibegin = temp.i - l1;
2682 iend = temp.i - 1;
2683 jbegin = temp.j;
2684 jend = temp.j + (int)strlen(temp.structure) - l1 - 2 - 1;
2685 for (k = ibegin + 1; k <= iend + 1; k++)
2686 struc2[k - 1] = temp.structure[k - ibegin - 1];
2687 for (k = jbegin + j; k <= jend + j; k++)
2688 struc2[k - 1] = temp.structure[l1 + k - j - jbegin + 1];
2689 }
2690
2691 free(temp.structure);
2692 }
2693
2694 strcpy(struc, st1 + MAX2(i - 1, 0));
2695 strcat(struc, "&");
2696 /* strcat(struc, st2); */
2697 strncat(struc, struc2 + 5, (int)strlen(struc2) - 10);
2698 free(struc2);
2699 free(struc_loop);
2700 free(st1);
2701 free(st2);
2702 /* free_arrays(); */
2703 return struc;
2704 }
2705
2706
2707 void
Lsnoop_subopt_list_XS(const char * s1,const char * s2,const int ** access_s1,int delta,int w,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int distance,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const int alignment_length,const char * name,const int fullStemEnergy)2708 Lsnoop_subopt_list_XS(const char *s1,
2709 const char *s2,
2710 const int **access_s1,
2711 int delta,
2712 int w,
2713 const int penalty,
2714 const int threshloop,
2715 const int threshLE,
2716 const int threshRE,
2717 const int threshDE,
2718 const int threshTE,
2719 const int threshSE,
2720 const int threshD,
2721 const int distance,
2722 const int half_stem,
2723 const int max_half_stem,
2724 const int min_s2,
2725 const int max_s2,
2726 const int min_s1,
2727 const int max_s1,
2728 const int min_d1,
2729 const int min_d2,
2730 const int alignment_length,
2731 const char *name,
2732 const int fullStemEnergy)
2733 {
2734 int min_colonne = INF;
2735 int max_pos;
2736 int max;
2737
2738 max = INF;
2739
2740 /* int temp; */
2741 /* int nsubopt=10; */
2742 n1 = (int)strlen(s1);
2743 n2 = (int)strlen(s2);
2744 int *position;
2745 int *position_j;
2746 int min_j_colonne;
2747 int max_pos_j = INF;
2748 position = (int *)vrna_alloc((n1 + 3) * sizeof(int));
2749 position_j = (int *)vrna_alloc((n1 + 3) * sizeof(int));
2750
2751 /* int Eminj, Emin_l; */
2752 int i, j;/* l1, Emin=INF, i_min=0, j_min=0; */
2753 /* char *struc; */
2754 /* snoopT mfe; */
2755 int *indx;
2756 int *mLoop;
2757 int *cLoop;
2758 folden **foldlist, **foldlist_XS;
2759 int Duplex_El, Duplex_Er;
2760 int Loop_D;
2761 /* int u; */
2762 int Loop_E;
2763 vrna_md_t md;
2764
2765 Duplex_El = 0;
2766 Duplex_Er = 0;
2767 Loop_E = 0, Loop_D = 0;
2768 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
2769
2770 set_model_details(&md);
2771 if ((!P) || (fabs(P->temperature - temperature) > 1e-6)) {
2772 snoupdate_fold_params();
2773 if (P)
2774 free(P);
2775
2776 P = vrna_params(&md);
2777 make_pair_matrix();
2778 }
2779
2780 lpair = (int **)vrna_alloc(sizeof(int *) * (6));
2781 lc = (int **)vrna_alloc(sizeof(int *) * (6));
2782 lr = (int **)vrna_alloc(sizeof(int *) * (6));
2783 for (i = 0; i < 6; i++) {
2784 lc[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
2785 lr[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
2786 lpair[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
2787 for (j = n2; j > -1; j--) {
2788 lc[i][j] = INF;
2789 lr[i][j] = INF;
2790 lpair[i][j] = 0;
2791 }
2792 }
2793 encode_seqs(s1, s2);
2794 int lim_maxj = n2 - min_d2;
2795 int lim_minj = min_d1;
2796 int lim_maxi = n1 - 5;
2797 for (i = 5; i <= lim_maxi; i++) {
2798 int idx = i % 5;
2799 int idx_1 = (i - 1) % 5;
2800 int idx_2 = (i - 2) % 5;
2801 int idx_3 = (i - 3) % 5;
2802 int idx_4 = (i - 4) % 5;
2803 int di1, di2, di3, di4;
2804 di1 = access_s1[5][i] - access_s1[4][i - 1];
2805 di2 = access_s1[5][i - 1] - access_s1[4][i - 2] + di1;
2806 di3 = access_s1[5][i - 2] - access_s1[4][i - 3] + di2;
2807 di4 = access_s1[5][i - 3] - access_s1[4][i - 4] + di3;
2808 di1 = MIN2(di1, 165);
2809 di2 = MIN2(di2, 330);
2810 di3 = MIN2(di3, 495);
2811 di4 = MIN2(di4, 660);
2812 for (j = lim_maxj; j > lim_minj; j--) {
2813 int type, type2;/* E, k,l; */
2814 type = pair[S1[i]][S2[j]];
2815 lpair[idx][j] = type;
2816 lc[idx][j] = (type) ? P->DuplexInit + access_s1[1][i] : INF;
2817 lr[idx][j] = INF;
2818 if (!type)
2819 continue;
2820
2821 if ( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
2822 j < max_s1 && j > min_s1 &&
2823 j > n2 - max_s2 - max_half_stem &&
2824 j < n2 - min_s2 - half_stem && S1[i - 2] == 4) {
2825 /*constraint on s2 and i*/
2826 int min_k, max_k;
2827 max_k = MIN2(n2 - min_s2, j + max_half_stem + 1);
2828 min_k = MAX2(j + half_stem + 1, n2 - max_s2);
2829 folden *temp;
2830 temp = foldlist[j + 1];
2831 while (temp->next) {
2832 int k = temp->k;
2833 /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
2834 if (lpair[idx_3][k + 1] && lc[idx_3][k + 1] /*+di3*/ < 411 /*&& lpair[idx_4][k+2]*/) /* remove second condition */
2835 lr[idx][j] = MIN2(lr[idx][j], di3 + lc[idx_3][k + 1] + temp->energy); /*--NU--*/
2836
2837 /*else*/ if (lpair[idx_4][k + 1] && /*di4 +*/ lc[idx_4][k + 1] < 411) /*--NUN--*/ /* remove second condition */
2838 lr[idx][j] = MIN2(lr[idx][j], di4 + lc[idx_4][k + 1] + temp->energy);
2839
2840 /* } */
2841 temp = temp->next;
2842 }
2843 }
2844
2845 /* dangle 5'SIDE relative to the mRNA */
2846 /**
2847 *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
2848 *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
2849 *** if (type>2) lc[idx][j] += P->TerminalAU;
2850 **/
2851 lc[idx][j] += vrna_E_ext_stem(type, SS1[i - 1], SS2[j + 1], P);
2852 /* if(j<n2 && i>1){ */
2853 /* type2=pair[S1[i-1]][S2[j+1]]; */
2854 type2 = lpair[idx_1][j + 1];
2855 if (type2 > 0) {
2856 lc[idx][j] =
2857 MIN2(lc[idx_1][j + 1] +
2858 E_IntLoop(0, 0, type2, rtype[type], SS1[i], SS2[j], SS1[i - 1], SS2[j + 1], P) + di1,
2859 lc[idx][j]);
2860 lr[idx][j] =
2861 MIN2(lr[idx_1][j + 1] +
2862 E_IntLoop(0, 0, type2, rtype[type], SS1[i], SS2[j], SS1[i - 1], SS2[j + 1], P) + di1,
2863 lr[idx][j]);
2864 }
2865
2866 type2 = lpair[idx_2][j + 2];
2867 if (type2 > 0) {
2868 lc[idx][j] =
2869 MIN2(lc[idx_2][j + 2] +
2870 E_IntLoop(1, 1, type2, rtype[type], SS1[i - 1], SS2[j + 1], SS1[i - 1], SS2[j + 1],
2871 P) + di2,
2872 lc[idx][j]);
2873 lr[idx][j] =
2874 MIN2(lr[idx_2][j + 2] +
2875 E_IntLoop(1, 1, type2, rtype[type], SS1[i - 1], SS2[j + 1], SS1[i - 1], SS2[j + 1],
2876 P) + di2,
2877 lr[idx][j]);
2878 }
2879
2880 type2 = lpair[idx_3][j + 3];
2881 if (type2 > 0) {
2882 lc[idx][j] =
2883 MIN2(lc[idx_3][j + 3] +
2884 E_IntLoop(2, 2, type2, rtype[type], SS1[i - 2], SS2[j + 2], SS1[i - 1], SS2[j + 1],
2885 P) + di3,
2886 lc[idx][j]);
2887 lr[idx][j] =
2888 MIN2(lr[idx_3][j + 3] +
2889 E_IntLoop(2, 2, type2, rtype[type], SS1[i - 2], SS2[j + 2], SS1[i - 1], SS2[j + 1],
2890 P) + di3,
2891 lr[idx][j]);
2892 }
2893
2894 int bla;
2895 int temp2;
2896 temp2 = min_colonne;
2897 bla = lr[idx][j] + vrna_E_ext_stem(rtype[type], SS2[j - 1], SS1[i + 1], P);
2898 /**
2899 *** (type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]];
2900 **/
2901 min_colonne = MIN2(bla, min_colonne);
2902 if (temp2 > min_colonne)
2903 min_j_colonne = j;
2904 }
2905 position[i] = min_colonne;
2906 if (max >= min_colonne) {
2907 max = min_colonne;
2908 max_pos = i;
2909 max_pos_j = min_j_colonne;
2910 }
2911
2912 position_j[i] = min_j_colonne;
2913 min_colonne = INF;
2914 }
2915 free(S1);
2916 free(S2);
2917 free(SS1);
2918 free(SS2);
2919
2920 if (max < threshTE + 30 * alignment_length) {
2921 find_max_snoop_XS(s1,
2922 s2,
2923 access_s1,
2924 max,
2925 alignment_length,
2926 position,
2927 position_j,
2928 delta,
2929 distance,
2930 penalty,
2931 threshloop,
2932 threshLE,
2933 threshRE,
2934 threshDE,
2935 threshTE,
2936 threshSE,
2937 threshD,
2938 half_stem,
2939 max_half_stem,
2940 min_s2,
2941 max_s2,
2942 min_s1,
2943 max_s1,
2944 min_d1,
2945 min_d2,
2946 name,
2947 fullStemEnergy);
2948 }
2949
2950 for (i = 1; i < 6; i++) {
2951 free(lc[i]);
2952 free(lr[i]);
2953 free(lpair[i]);
2954 }
2955 free(lc[0]);
2956 free(lr[0]);
2957 free(lpair[0]);
2958 free(lc);
2959 free(lr);
2960 free(lpair);
2961 free(position);
2962 free(position_j);
2963 }
2964
2965
2966 PRIVATE void
find_max_snoop_XS(const char * s1,const char * s2,const int ** access_s1,const int max,const int alignment_length,const int * position,const int * position_j,const int delta,const int distance,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshTE,const int threshSE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const char * name,const int fullStemEnergy)2967 find_max_snoop_XS(const char *s1,
2968 const char *s2,
2969 const int **access_s1,
2970 const int max,
2971 const int alignment_length,
2972 const int *position,
2973 const int *position_j,
2974 const int delta,
2975 const int distance,
2976 const int penalty,
2977 const int threshloop,
2978 const int threshLE,
2979 const int threshRE,
2980 const int threshDE,
2981 const int threshTE,
2982 const int threshSE,
2983 const int threshD,
2984 const int half_stem,
2985 const int max_half_stem,
2986 const int min_s2,
2987 const int max_s2,
2988 const int min_s1,
2989 const int max_s1,
2990 const int min_d1,
2991 const int min_d2,
2992 const char *name,
2993 const int fullStemEnergy)
2994 {
2995 int count = 0;
2996 int n3 = (int)strlen(s1);
2997 int n4 = (int)strlen(s2);
2998 int pos = n1 - 4;
2999 int max_pos_j;
3000 int threshold = MIN2(threshTE + alignment_length * 30, -100);
3001
3002 /* printf("threshTE %d max %d\n", threshTE, max); */
3003 /* #pragma omp parallel for */
3004 /* for(pos=n1+1;pos>distance;pos--){ */
3005 while (pos-- > 5) {
3006 int temp_min = 0;
3007 if (position[pos] < (threshold)) {
3008 int search_range;
3009 search_range = distance + 1;
3010 while (--search_range)
3011 if (position[pos - search_range] <= position[pos - temp_min])
3012 temp_min = search_range;
3013
3014 pos -= temp_min;
3015 max_pos_j = position_j[pos];
3016 int begin = MAX2(5, pos - alignment_length);
3017 int end = MIN2(n3 - 5, pos - 1);
3018 char *s3 = (char *)vrna_alloc(sizeof(char) * (end - begin + 2) + 5);
3019 strncpy(s3, (s1 + begin), end - begin + 1);
3020 strcat(s3, "NNNNN\0");
3021
3022 int n5 = (int)strlen(s3);
3023 snoopT test;
3024 test = snoopfold_XS(s3, s2, access_s1, pos, max_pos_j, penalty,
3025 threshloop, threshLE, threshRE,
3026 threshDE, threshD, half_stem,
3027 max_half_stem, min_s2, max_s2, min_s1,
3028 max_s1, min_d1, min_d2, fullStemEnergy);
3029 if (test.energy == INF) {
3030 free(s3);
3031 continue;
3032 }
3033
3034 if (test.Duplex_El > threshLE * 0.01 || test.Duplex_Er > threshRE * 0.01 ||
3035 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
3036 (test.Duplex_Er + test.Duplex_El + test.Loop_E) > threshTE * 0.01 ||
3037 (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE * 0.01) {
3038 free(test.structure);
3039 free(s3);
3040 continue;
3041 }
3042
3043 char *s4;
3044 s4 = (char *)vrna_alloc(sizeof(char) * (n4 - 9));
3045 strncpy(s4, s2 + 5, n4 - 10);
3046 s4[n4 - 10] = '\0';
3047
3048 char *s5 = vrna_alloc(sizeof(char) * n5 - test.i + 2 - 5);
3049 strncpy(s5, s3 + test.i - 1, n5 - test.i + 1 - 5);
3050 s5[n5 - test.i + 1 - 5] = '\0';
3051 float dE = ((float)(access_s1[n5 - test.i + 1 - 5][pos])) * 0.01;
3052 printf(
3053 "%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10) (%5.2f)\n%s&%s\n",
3054 test.structure,
3055 pos - (n5 - test.i),
3056 pos - 5,
3057 pos - (n5 - test.u),
3058 max_pos_j - 5,
3059 max_pos_j - 5 + (int)(strrchr(test.structure, '>') - strchr(test.structure, '>')),
3060 test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10 + dE,
3061 test.Duplex_El,
3062 test.Duplex_Er,
3063 test.Loop_E,
3064 test.Loop_D,
3065 dE,
3066 test.fullStemEnergy,
3067 s5,
3068 s4);
3069 if (name) {
3070 int begin_t, end_t, begin_q, end_q, and, pipe, i;
3071 char *psoutput;
3072 begin_q = 0;
3073 end_q = n4 - 10;
3074 begin_t = 0;
3075 end_t = n5 - test.i + 1 - 5;
3076 and = end_t + 1;
3077 pipe = test.u - test.i + 1;
3078 cut_point = end_t + 1;
3079 char *catseq, *catstruct;/* *fname; */
3080 catseq = (char *)vrna_alloc(n5 + end_q - begin_q + 2);
3081 catstruct = (char *)vrna_alloc(n5 + end_q - begin_q + 2);
3082 strcpy(catseq, s5);
3083 strncpy(catstruct, test.structure, end_t);
3084 strcat(catseq, s4);
3085 strncat(catstruct, test.structure + end_t + 1, end_q - begin_q + 1);
3086 catstruct[end_t - begin_t + end_q - begin_q + 2] = '\0';
3087 catseq[end_t - begin_t + end_q - begin_q + 2] = '\0';
3088 int *relative_access;
3089 relative_access = vrna_alloc(sizeof(int) * strlen(s5));
3090
3091 relative_access[0] = access_s1[1][pos - (n5 - test.i) + 5];
3092 for (i = 1; i < (int)strlen(s5); i++)
3093 relative_access[i] = access_s1[i + 1][pos - (n5 - test.i) + i + 5] -
3094 access_s1[i][pos - (n5 - test.i) + i + 4];
3095
3096 psoutput = vrna_strdup_printf("sno_XS_%d_u_%d_%s.ps",
3097 count,
3098 pos - (n5 - test.u),
3099 name);
3100
3101 PS_rna_plot_snoop_a(catseq, catstruct, psoutput, relative_access, NULL);
3102 free(catseq);
3103 free(catstruct);
3104 free(relative_access);
3105 free(psoutput);
3106 count++;
3107 }
3108
3109 free(s3);
3110 free(s4);
3111 free(s5);
3112 free(test.structure);
3113 }
3114 }
3115 }
3116
3117
3118 snoopT
snoopfold_XS(const char * s1,const char * s2,const int ** access_s1,const int pos_i,const int pos_j,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2,const int fullStemEnergy)3119 snoopfold_XS(const char *s1,
3120 const char *s2,
3121 const int **access_s1,
3122 const int pos_i,
3123 const int pos_j,
3124 const int penalty,
3125 const int threshloop,
3126 const int threshLE,
3127 const int threshRE,
3128 const int threshDE,
3129 const int threshD,
3130 const int half_stem,
3131 const int max_half_stem,
3132 const int min_s2,
3133 const int max_s2,
3134 const int min_s1,
3135 const int max_s1,
3136 const int min_d1,
3137 const int min_d2,
3138 const int fullStemEnergy)
3139 {
3140 /* int Eminj, Emin_l; */
3141 int a, b, i, j, Emin = INF, a_min = 0, b_min = 0;
3142 char *struc;
3143 snoopT mfe;
3144 int *indx;
3145 int *mLoop;
3146 int *cLoop;
3147 folden **foldlist, **foldlist_XS;
3148 int Duplex_El, Duplex_Er;
3149 int Loop_D;
3150 int u;
3151 int Loop_E;
3152 vrna_md_t md;
3153
3154 Duplex_El = 0;
3155 Duplex_Er = 0;
3156 Loop_E = 0, Loop_D = 0;
3157 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
3158 n1 = (int)strlen(s1);
3159 n2 = (int)strlen(s2);
3160
3161 set_model_details(&md);
3162 if ((!P) || (fabs(P->temperature - temperature) > 1e-6)) {
3163 snoupdate_fold_params();
3164 if (P)
3165 free(P);
3166
3167 P = vrna_params(&md);
3168 make_pair_matrix();
3169 }
3170
3171 c = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
3172 r = (int **)vrna_alloc(sizeof(int *) * (n1 + 1));
3173 for (i = 0; i <= n1; i++) {
3174 c[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
3175 r[i] = (int *)vrna_alloc(sizeof(int) * (n2 + 1));
3176 for (j = n2; j > -1; j--) {
3177 c[i][j] = INF;
3178 r[i][j] = INF;
3179 }
3180 }
3181 encode_seqs(s1, s2);
3182 i = n1 - 5;
3183 j = pos_j;
3184 /* printf("tar: %s\nsno: %s\n ", s1, s2); */
3185 /* printf("pos_i %d pos_j %d\n", pos_i, pos_j); */
3186 /* printf("type %d n1 %d n2 %d S1[n1] %d S2[n2] %d", pair[S1[i]][S2[j]], n1, n2, S1[n1], S2[n2]); */
3187 int type, type2, E, p, q;
3188 r[i][j] = P->DuplexInit;
3189 /* r[i][j] += P->dangle3[rtype[type]][SS1[i+1]] + P->dangle5[rtype[type]][SS2[j-1]]; */
3190
3191 if (pair[S1[i]][S2[j]] > 2)
3192 r[i][j] += P->TerminalAU;
3193
3194 for (a = i - 1; a > 0; a--) {
3195 /* i-1 */
3196 r[a + 1][0] = INF;
3197 for (b = j + 1; b <= n2 - min_d2; b++) {
3198 /* j+1 */
3199 r[a][b] = INF;
3200 type = pair[S1[a]][S2[b]];
3201 if (!type)
3202 continue;
3203
3204 if (S1[a + 1] == 4) {
3205 folden *temp;
3206 temp = foldlist_XS[b - 1];
3207 while (temp->next) {
3208 int k = temp->k;
3209 if (pair[S1[a + 3]][S2[k - 1]] && k < max_s1 && k > min_s1 &&
3210 k > n2 - max_s2 - max_half_stem &&
3211 k < n2 - min_s2 - half_stem /*&& r[a+3][k-1] + access_s1[i-(a+3)+1][pos_i] < 411*/) /* remove last condition last condition is to check if the interaction is stable enough */
3212 c[a][b] = MIN2(c[a][b], r[a + 3][k - 1] + temp->energy);
3213
3214 temp = temp->next;
3215 }
3216 }
3217
3218 if (S1[a + 2] == 4) {
3219 folden *temp;
3220 temp = foldlist_XS[b - 1];
3221 while (temp->next) {
3222 int k = temp->k;
3223 if (pair[S1[a + 4]][S2[k - 1]] && k < max_s1 && k > min_s1 &&
3224 k > n2 - max_s2 - max_half_stem &&
3225 k < n2 - min_s2 - half_stem /*&& r[a+4][k-1] + access_s1[i-(a+4)+1][pos_i] < 411 */) /* remove last condition */
3226 c[a][b] = MIN2(c[a][b], r[a + 4][k - 1] + temp->energy);
3227
3228 temp = temp->next;
3229 }
3230 }
3231
3232 for (p = a + 1; p < n1 && (p - a) < MAXLOOP_L; p++) {
3233 /* p < n1 */
3234 for (q = b - 1; q > 1; q--) {
3235 /* q > 1 */
3236 if (p - a + b - q > 2 * MAXLOOP_L - 2)
3237 break;
3238
3239 if (abs((p - a) - (b - q)) >= ASS)
3240 continue;
3241
3242 type2 = pair[S1[p]][S2[q]];
3243 if (!type2)
3244 continue;
3245
3246 E =
3247 E_IntLoop(p - a - 1,
3248 b - q - 1,
3249 type2,
3250 rtype[type],
3251 SS1[a + 1],
3252 SS2[b - 1],
3253 SS1[p - 1],
3254 SS2[q + 1],
3255 P);
3256 c[a][b] = MIN2(c[a][b], c[p][q] + E);
3257 r[a][b] = MIN2(r[a][b], r[p][q] + E);
3258 }
3259 }
3260 E = c[a][b];
3261 if (type > 2)
3262 E += P->TerminalAU;
3263
3264 /* E +=P->dangle5[rtype[type]][SS1[i+1]]; */
3265 /* E +=P->dangle5[rtype[type]][SS2[j-1]]; */
3266 E += access_s1[i - a + 1][pos_i];
3267 if (E < Emin) {
3268 Emin = E;
3269 a_min = a;
3270 b_min = b;
3271 }
3272 }
3273 }
3274 if (Emin > 0) {
3275 printf("no target found under the constraints chosen\n");
3276 for (i = 0; i <= n1; i++) {
3277 free(r[i]);
3278 free(c[i]);
3279 }
3280 free(c);
3281 free(r);
3282 free(S1);
3283 free(S2);
3284 free(SS1);
3285 free(SS2);
3286 mfe.energy = INF;
3287 return mfe;
3288 }
3289
3290 type2 = pair[S1[a_min]][S2[b_min]];
3291 if (type2 > 2)
3292 Emin += P->TerminalAU;
3293
3294 mfe.energy = ((float)(Emin)) / 100;
3295 struc = snoop_backtrack_XS(a_min,
3296 b_min,
3297 s2,
3298 &Duplex_El,
3299 &Duplex_Er,
3300 &Loop_E,
3301 &Loop_D,
3302 &u,
3303 penalty,
3304 threshloop,
3305 threshLE,
3306 threshRE,
3307 threshDE,
3308 threshD,
3309 half_stem,
3310 max_half_stem,
3311 min_s2,
3312 max_s2,
3313 min_s1,
3314 max_s1,
3315 min_d1,
3316 min_d2);
3317
3318 mfe.i = a_min;
3319 mfe.j = b_min;
3320 mfe.u = u;
3321 mfe.Duplex_Er = (float)Duplex_Er / 100;
3322 mfe.Duplex_El = (float)Duplex_El / 100;
3323 mfe.Loop_D = (float)Loop_D / 100;
3324 mfe.Loop_E = (float)Loop_E / 100;
3325 mfe.energy = (float)Emin / 100;
3326 mfe.fullStemEnergy = (float)fullStemEnergy / 100;
3327 mfe.structure = struc;
3328 return mfe;
3329 }
3330
3331
3332 PRIVATE char *
snoop_backtrack_XS(int i,int j,const char * snoseq,int * Duplex_El,int * Duplex_Er,int * Loop_E,int * Loop_D,int * u,const int penalty,const int threshloop,const int threshLE,const int threshRE,const int threshDE,const int threshD,const int half_stem,const int max_half_stem,const int min_s2,const int max_s2,const int min_s1,const int max_s1,const int min_d1,const int min_d2)3333 snoop_backtrack_XS(int i,
3334 int j,
3335 const char *snoseq,
3336 int *Duplex_El,
3337 int *Duplex_Er,
3338 int *Loop_E,
3339 int *Loop_D,
3340 int *u,
3341 const int penalty,
3342 const int threshloop,
3343 const int threshLE,
3344 const int threshRE,
3345 const int threshDE,
3346 const int threshD,
3347 const int half_stem,
3348 const int max_half_stem,
3349 const int min_s2,
3350 const int max_s2,
3351 const int min_s1,
3352 const int max_s1,
3353 const int min_d1,
3354 const int min_d2)
3355 {
3356 /* backtrack structure going backwards from i, and forwards from j
3357 * return structure in bracket notation with & as separator */
3358 int k, l, type, type2, E, traced, i0, j0;
3359 int traced_c = 0; /* flag for following backtrack in c or r */
3360 char *st1, *st2, *struc;
3361 char *struc_loop;
3362
3363 st1 = (char *)vrna_alloc(sizeof(char) * (n1 + 1));
3364 st2 = (char *)vrna_alloc(sizeof(char) * (n2 + 1));
3365 int *indx;
3366 int *mLoop;
3367 int *cLoop;
3368 folden **foldlist, **foldlist_XS;
3369 type = pair[S1[i]][S2[j]];
3370 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
3371 i0 = i;
3372 j0 = j;
3373 /* i0=MAX2(i,1); j0=MIN2(j+1,n2); */
3374 while (i <= n1 && j >= 1) {
3375 if (!traced_c) {
3376 E = c[i][j];
3377 traced = 0;
3378 st1[i] = '<';
3379 st2[j - 1] = '>';
3380 type = pair[S1[i]][S2[j]];
3381 if (!type)
3382 vrna_message_error("backtrack failed in fold duplex c");
3383
3384 for (k = i + 1; k > 0 && (k - i) < MAXLOOP_L; k++) {
3385 for (l = j - 1; l >= 1; l--) {
3386 int LE;
3387 if (k - i + j - l > 2 * MAXLOOP_L - 2)
3388 break;
3389
3390 if (abs(k - i - j + l) >= ASS)
3391 continue;
3392
3393 type2 = pair[S1[k]][S2[l]];
3394 if (!type2)
3395 continue;
3396
3397 LE = E_IntLoop(k - i - 1, j - l - 1, type2, rtype[type],
3398 SS1[i + 1], SS2[j - 1], SS1[k - 1], SS2[l + 1], P);
3399 if (E == c[k][l] + LE) {
3400 traced = 1;
3401 i = k;
3402 j = l;
3403 *Duplex_El += LE;
3404 break;
3405 }
3406 }
3407 if (traced)
3408 break;
3409 }
3410 if (!traced) {
3411 if (S1[i + 1] == 4) {
3412 folden *temp;
3413 temp = foldlist_XS[j - 1];
3414 while (temp->next) {
3415 int k = temp->k;
3416 if (pair[S1[i + 3]][S2[k - 1]] && k < max_s1 && k > min_s1 &&
3417 k > n2 - max_s2 - max_half_stem && k < n2 - min_s2 - half_stem) {
3418 if (E == r[i + 3][k - 1] + temp->energy) {
3419 *Loop_E = temp->energy;
3420 st1[i + 1] = '|';
3421 st1[i + 2] = '.';
3422 *u = i + 1;
3423 int a, b;
3424 for (a = 0; a < MISMATCH; a++) {
3425 for (b = 0; b < MISMATCH; b++) {
3426 int ij = indx[j - 1 - a] + k + b;
3427 if (cLoop[ij] == temp->energy) {
3428 struc_loop = snobacktrack_fold_from_pair(snoseq, k + b, j - 1 - a);
3429 a = INF;
3430 b = INF;
3431 }
3432 }
3433 }
3434 traced = 1;
3435 traced_c = 1;
3436 i = i + 3;
3437 j = k - 1;
3438 break;
3439 }
3440 }
3441
3442 temp = temp->next;
3443 }
3444 }
3445
3446 if (S1[i + 2] == 4) {
3447 /* introduce structure from RNAfold */
3448 folden *temp;
3449 temp = foldlist_XS[j - 1];
3450 while (temp->next) {
3451 int k = temp->k;
3452 if (pair[S1[i + 4]][S2[k - 1]] && k < max_s1 && k > min_s1 &&
3453 k > n2 - max_s2 - max_half_stem && k < n2 - min_s2 - half_stem) {
3454 if (E == r[i + 4][k - 1] + temp->energy) {
3455 *Loop_E = temp->energy;
3456 st1[i + 2] = '|';
3457 st1[i + 1] = st1[i + 3] = '.';
3458 *u = i + 2;
3459 int a, b;
3460 for (a = 0; a < MISMATCH; a++) {
3461 for (b = 0; b < MISMATCH; b++) {
3462 int ij = indx[j - 1 - a] + k + b;
3463 if (cLoop[ij] == temp->energy) {
3464 struc_loop = snobacktrack_fold_from_pair(snoseq, k + b, j - a - 1);
3465 a = INF;
3466 b = INF;
3467 }
3468 }
3469 }
3470 traced = 1;
3471 traced_c = 1;
3472 i = i + 4;
3473 j = k - 1;
3474 break;
3475 }
3476 }
3477
3478 temp = temp->next;
3479 }
3480 }
3481 } /* traced? */
3482 } /* traced_r? */
3483 else {
3484 E = r[i][j];
3485 traced = 0;
3486 st1[i] = '<';
3487 st2[j - 1] = '>';
3488 type = pair[S1[i]][S2[j]];
3489 if (!type)
3490 vrna_message_error("backtrack failed in fold duplex r");
3491
3492 for (k = i + 1; k > 0 && (k - i) < MAXLOOP_L; k++) {
3493 for (l = j - 1; l >= 1; l--) {
3494 int LE;
3495 if (k - i + j - l > 2 * MAXLOOP_L - 2)
3496 break;
3497
3498 if (abs(k - i - j + l) >= ASS)
3499 continue;
3500
3501 type2 = pair[S1[k]][S2[l]];
3502 if (!type2)
3503 continue;
3504
3505 LE = E_IntLoop(k - i - 1, j - l - 1, type2, rtype[type],
3506 SS1[i + 1], SS2[j - 1], SS1[k - 1], SS2[l + 1], P);
3507 if (E == r[k][l] + LE) {
3508 traced = 1;
3509 i = k;
3510 j = l;
3511 *Duplex_Er += LE;
3512 break;
3513 }
3514 }
3515 if (traced)
3516 break;
3517 }
3518 }
3519
3520 if (!traced) {
3521 /* if (i>1) {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];} */
3522 /* if (j<n2) {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];} */
3523 if (type > 2) {
3524 E -= P->TerminalAU;
3525 *Duplex_Er += P->TerminalAU;
3526 }
3527
3528 if (E != P->DuplexInit)
3529 vrna_message_error("backtrack failed in fold duplex end");
3530 else
3531 break;
3532 }
3533 }
3534
3535
3536 /* struc = (char *) vrna_alloc(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
3537 struc = (char *)vrna_alloc(i - i0 + 1 + n2); /* declare final duplex structure */
3538 char *struc2;
3539 struc2 = (char *)vrna_alloc(n2 + 1);
3540 /* char * struct_const; */
3541
3542 for (k = MIN2(i0, 1); k <= i; k++)
3543 if (!st1[k - 1])
3544 st1[k - 1] = '.';
3545
3546 /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
3547 /* char * struct_const; */
3548 /* struct_const = (char *) vrna_alloc(sizeof(char)*(n2+1)); */
3549 for (k = 1; k <= n2; k++) {
3550 if (!st2[k - 1])
3551 st2[k - 1] = struc_loop[k - 1]; /* '.'; */
3552
3553 struc2[k - 1] = st2[k - 1]; /* '.'; */
3554 /* if (k>=j0 && k<=j){ */
3555 /* struct_const[k-1]='x'; */
3556 /* } */
3557 /* else{ */
3558 /* if(k<j0) {struct_const[k-1]='<';} */
3559 /* if(k>j) {struct_const[k-1]='>';} */
3560 /* } */
3561 }
3562 char duplexseq_1[j];
3563 char duplexseq_2[n2 - j0 + 2];
3564 if (j0 < n2) {
3565 strncpy(duplexseq_1, snoseq, j - 1);
3566 strcpy(duplexseq_2, snoseq + j0);
3567 duplexseq_1[j - 1] = '\0';
3568 duplexseq_2[n2 - j0 + 1] = '\0';
3569 duplexT temp;
3570 temp = duplexfold(duplexseq_1, duplexseq_2);
3571 *Loop_D = MIN2(0, -410 + (int)100 * temp.energy);
3572 if (*Loop_D) {
3573 int l1, ibegin, iend, jbegin, jend;
3574 l1 = strchr(temp.structure, '&') - temp.structure;
3575 ibegin = temp.i - l1;
3576 iend = temp.i - 1;
3577 jbegin = temp.j;
3578 jend = temp.j + (int)strlen(temp.structure) - l1 - 2 - 1;
3579 for (k = ibegin + 1; k <= iend + 1; k++)
3580 struc2[k - 1] = temp.structure[k - ibegin - 1];
3581 for (k = jbegin + j0; k <= jend + j0; k++)
3582 struc2[k - 1] = temp.structure[l1 + k - j0 - jbegin + 1];
3583 }
3584
3585 free(temp.structure);
3586 }
3587
3588 strcpy(struc, st1 + MAX2(i0, 1));
3589 strcat(struc, "&");
3590 /* strcat(struc, st2); */
3591 strncat(struc, struc2 + 5, (int)strlen(struc2) - 10);
3592 free(struc2);
3593 free(struc_loop);
3594 free(st1);
3595 free(st2);
3596
3597 for (i = 0; i <= n1; i++) {
3598 free(r[i]);
3599 free(c[i]);
3600 }
3601 free(c);
3602 free(r);
3603 free(S1);
3604 free(S2);
3605 free(SS1);
3606 free(SS2);
3607 /* free_arrays(); */
3608 return struc;
3609 }
3610
3611
3612 PRIVATE int
covscore(const int * types,int n_seq)3613 covscore(const int *types,
3614 int n_seq)
3615 {
3616 /* calculate co-variance bonus for a pair depending on */
3617 /* compensatory/consistent mutations and incompatible seqs */
3618 /* should be 0 for conserved pairs, >0 for good pairs */
3619 #define NONE -10000 /* score for forbidden pairs */
3620 int k, l, s, score, pscore;
3621 int dm[7][7] = { { 0, 0, 0, 0, 0, 0, 0 }, /* hamming distance between pairs */
3622 { 0, 0, 2, 2, 1, 2, 2 } /* CG */,
3623 { 0, 2, 0, 1, 2, 2, 2 } /* GC */,
3624 { 0, 2, 1, 0, 2, 1, 2 } /* GU */,
3625 { 0, 1, 2, 2, 0, 2, 1 } /* UG */,
3626 { 0, 2, 2, 1, 2, 0, 2 } /* AU */,
3627 { 0, 2, 2, 2, 1, 2, 0 } /* UA */ };
3628
3629 int pfreq[8] = {
3630 0, 0, 0, 0, 0, 0, 0, 0
3631 };
3632 for (s = 0; s < n_seq; s++)
3633 pfreq[types[s]]++;
3634
3635 if (pfreq[0] * 2 > n_seq)
3636 return NONE;
3637
3638 for (k = 1, score = 0; k <= 6; k++) /* ignore pairtype 7 (gap-gap) */
3639 for (l = k + 1; l <= 6; l++)
3640 /* scores for replacements between pairtypes */
3641 /* consistent or compensatory mutations score 1 or 2 */
3642 score += pfreq[k] * pfreq[l] * dm[k][l];
3643
3644 /* counter examples score -1, gap-gap scores -0.25 */
3645 pscore = cv_fact *
3646 ((UNIT * score) / n_seq - nc_fact * UNIT * (pfreq[0] + pfreq[7] * 0.25));
3647 return pscore;
3648 }
3649
3650
3651 /*---------------------------------------------------------------------------*/
3652
3653 PRIVATE short *
aliencode_seq(const char * sequence)3654 aliencode_seq(const char *sequence)
3655 {
3656 unsigned int i, l;
3657 short *Stemp;
3658
3659 l = strlen(sequence);
3660 Stemp = (short *)vrna_alloc(sizeof(short) * (l + 2));
3661 Stemp[0] = (short)l;
3662
3663 /* make numerical encoding of sequence */
3664 for (i = 1; i <= l; i++)
3665 Stemp[i] = (short)encode_char(toupper(sequence[i - 1]));
3666
3667 /* for circular folding add first base at position n+1 */
3668 /* Stemp[l+1] = Stemp[1]; */
3669
3670 return Stemp;
3671 }
3672
3673
3674 PRIVATE short *
encode_seq(const char * sequence)3675 encode_seq(const char *sequence)
3676 {
3677 unsigned int i, l;
3678 short *S;
3679
3680 l = strlen(sequence);
3681 extern double nc_fact;
3682 S = (short *)vrna_alloc(sizeof(short) * (l + 2));
3683 S[0] = (short)l;
3684
3685 /* make numerical encoding of sequence */
3686 for (i = 1; i <= l; i++)
3687 S[i] = (short)encode_char(toupper(sequence[i - 1]));
3688 /* for circular folding add first base at position n+1 */
3689 S[l + 1] = S[1];
3690
3691 return S;
3692 }
3693
3694
3695 PRIVATE void
encode_seqs(const char * s1,const char * s2)3696 encode_seqs(const char *s1,
3697 const char *s2)
3698 {
3699 unsigned int i, l;
3700
3701 l = strlen(s1);
3702 S1 = encode_seq(s1);
3703 SS1 = (short *)vrna_alloc(sizeof(short) * (l + 1));
3704 /* SS1 exists only for the special X K and I bases and energy_set!=0 */
3705
3706 for (i = 1; i <= l; i++) /* make numerical encoding of sequence */
3707 SS1[i] = alias[S1[i]]; /* for mismatches of nostandard bases */
3708
3709 l = strlen(s2);
3710 S2 = encode_seq(s2);
3711 SS2 = (short *)vrna_alloc(sizeof(short) * (l + 1));
3712 /* SS2 exists only for the special X K and I bases and energy_set!=0 */
3713
3714 for (i = 1; i <= l; i++) /* make numerical encoding of sequence */
3715 SS2[i] = alias[S2[i]]; /* for mismatches of nostandard bases */
3716 }
3717
3718
3719 PRIVATE int
compare(const void * sub1,const void * sub2)3720 compare(const void *sub1,
3721 const void *sub2)
3722 {
3723 int d;
3724
3725 if (((snoopT *)sub1)->energy > ((snoopT *)sub2)->energy)
3726 return 1;
3727
3728 if (((snoopT *)sub1)->energy < ((snoopT *)sub2)->energy)
3729 return -1;
3730
3731 d = ((snoopT *)sub1)->i - ((snoopT *)sub2)->i;
3732 if (d != 0)
3733 return d;
3734
3735 return ((snoopT *)sub1)->j - ((snoopT *)sub2)->j;
3736 }
3737