1 /* Last changed Time-stamp: <2008-07-04 15:57:03 ulim> */
3 partiton function for RNA secondary structures
10 $Log: part_func_up.c,v $
11 Revision 1.4 2008/07/04 14:27:36 ivo
14 Revision 1.3 2008/05/08 14:11:55 ivo
17 Revision 1.2 2007/12/13 10:19:54 ivo
18 major RNAup update from Ulli
20 Revision 1.1 2007/04/30 15:13:13 ivo
21 merge RNAup into package
23 Revision 1.11 2006/07/17 11:11:43 ulim
24 removed all globals from fold_vars.h,c, cleaned code
26 Revision 1.10 2006/07/12 09:19:29 ulim
27 global variables w, incr3 and incr5 are now local
29 Revision 1.9 2006/07/11 12:45:02 ulim
30 remove redundancy in function pf_interact(...)
32 Revision 1.8 2006/03/08 15:26:37 ulim
33 modified -o[1|2], added meaningful default
35 Revision 1.5 2006/01/23 11:27:04 ulim
36 include file into new package version. cleaned it
38 Revision 1.2 2005/07/29 15:13:37 ulim
39 put the function, calculating the probability of an unpaired region in
40 an RNA and the function calculating the prob. of interaction between 2 RNAs
41 in a seperate file (pf_two.c)
43 Revision 1.1 2005/07/26 13:27:12 ulim
46 Revision 1.2 2005/07/01 13:14:57 ulim
47 fixed error in scaling, included new commandline options -incr5, -incr3 to
48 allow a variable number of unpaired positions 5' and 3' of the site of
49 interaction between the two RNAs
51 Revision 1.1 2005/04/19 08:16:38 ulim
60 #include <float.h> /* #defines FLT_MAX ... */
64 #include "energy_par.h"
65 #include "fold_vars.h"
68 #include "part_func.h"
69 #include "loop_energies.h"
70 #include "part_func_up.h"
73 static char rcsid[] UNUSED = "$Id: part_func_up.c,v 1.4 2008/07/04 14:27:36 ivo Exp $";
76 #define ZERO(A) (fabs(A) < DBL_EPSILON)
77 #define EQUAL(A,B) (fabs((A)-(B)) < 1000*DBL_EPSILON)
78 #define ISOLATED 256.0
79 /* #define NUMERIC 1 */
82 #################################
84 #################################
88 #################################
90 #################################
92 PRIVATE short *S=NULL, *S1=NULL, *SS=NULL, *SS2=NULL;
93 PRIVATE pf_paramT *Pf = NULL;/* use this structure for all the exp-arrays*/
94 PRIVATE FLT_OR_DBL *qb=NULL, *qm=NULL, *prpr=NULL; /* add arrays for pf_unpaired()*/
95 PRIVATE FLT_OR_DBL *probs=NULL;
96 PRIVATE FLT_OR_DBL *q1k=NULL, *qln=NULL;
97 PRIVATE double *qqm2=NULL, *qq_1m2=NULL, *qqm=NULL, *qqm1=NULL;
98 PRIVATE FLT_OR_DBL *scale=NULL, *expMLbase=NULL;
99 PRIVATE char *ptype=NULL; /* precomputed array of pair types */
100 PRIVATE int init_length; /* length in last call to init_pf_fold()*/
101 PRIVATE double init_temp; /* temperature in last call to scale_pf_params */
102 PRIVATE int *my_iindx = NULL;
103 /* make iptypes array for intermolecular constrains (ipidx for indexing)*/
107 #################################
108 # PRIVATE FUNCTION DECLARATIONS #
109 #################################
111 PRIVATE pu_out *get_u_vals(pu_contrib *p_c,
112 int **unpaired_values,
113 char *select_contrib);
115 PRIVATE int plot_free_pu_out( pu_out* res,
120 PRIVATE void scale_stru_pf_params(unsigned int length);
122 PRIVATE void init_pf_two(int length);
124 PRIVATE void scale_int(const char *s,
128 PRIVATE void encode_seq( const char *s1,
131 PRIVATE constrain *get_ptypes(char *S,
132 const char *structure);
134 PRIVATE void get_up_arrays(unsigned int length);
136 PRIVATE void free_up_arrays(void);
138 PRIVATE void set_encoded_seq(const char *sequence,
142 PRIVATE void get_interact_arrays(unsigned int n1,
153 #################################
154 # BEGIN OF FUNCTION DEFINITIONS #
155 #################################
158 PUBLIC pu_contrib *get_pu_contrib_struct(unsigned int n, unsigned int w){
160 pu_contrib *pu = (pu_contrib *)space(sizeof(pu_contrib));
163 /* contributions to probability of being unpaired witihin a(n)
168 /* pu_test->X[i][j] where i <= j and i [1...n], j = [1...w[ */
169 pu->H = (double **)space(sizeof(double *) * (n + 1));
170 pu->I = (double **)space(sizeof(double *) * (n + 1));
171 pu->M = (double **)space(sizeof(double *) * (n + 1));
172 pu->E = (double **)space(sizeof(double *) * (n + 1));
174 pu->H[i] = (double *)space(sizeof(double) * (w + 1));
175 pu->I[i] = (double *)space(sizeof(double) * (w + 1));
176 pu->M[i] = (double *)space(sizeof(double) * (w + 1));
177 pu->E[i] = (double *)space(sizeof(double) * (w + 1));
182 PUBLIC void free_pu_contrib_struct(pu_contrib *pu){
185 for(i=0;i<=pu->length;i++){
199 /* you have to call pf_fold(sequence, structure); befor pf_unstru */
200 PUBLIC pu_contrib *pf_unstru(char *sequence, int w){
201 int n, i, j, v, k, l, o, p, ij, kl, po, u, u1, d, type, type_2, tt, uu;
204 double qbt1, *tmp, Zuij, sum_l, *sum_M;
205 double *store_H, *store_Io, **store_I2o; /* hairp., interior contribs */
206 double *store_M_qm_o,*store_M_mlbase; /* multiloop contributions */
211 n = (int) strlen(sequence);
212 sum_M = (double *) space((n+1) * sizeof(double));
213 pu_test = get_pu_contrib_struct((unsigned)n, (unsigned)w);
214 size = ((n+1)*(n+2))>>1;
216 get_up_arrays((unsigned) n);
219 /* init everything */
220 for (d=0; d<=TURN; d++)
221 for (i=1; i<=n-d; i++){
225 pu_test->H[i][d]=pu_test->I[i][d]=pu_test->M[i][d]=pu_test->E[i][d]=0.;
230 for (i=0; i<size; i++)
234 for (i=1; i<=n; i++){
235 /* set auxillary arrays to 0, reuse qqm and qqm1, reuse qqm2 and qq_1m2*/
236 sum_M[i] = qqm[i] = qqm1[i] = qqm2[i] = qq_1m2[i] = 0;
237 for (j=i+TURN+1; j<=n; j++){
239 /* i need the part_func of all structures outside bp[ij] */
240 if(qb[ij] > 0.0) prpr[ij]= (probs[ij]/qb[ij]);
244 /* alloc even more memory */
245 store_I2o = (double **)space(sizeof(double *) * (n + 1)); /* for p,k */
247 store_I2o[i] = (double *)space(sizeof(double) * (MAXLOOP + 2));
249 /* expMLbase[i-p]*dangles_po */
250 store_M_mlbase = (double *)space(sizeof(double) * (size + 1));
252 /* 2. exterior bp (p,o) encloses unpaired region [i,i+w[*/
253 for (o=TURN+2;o<=n; o++) {
255 /*allocate space for arrays to store different contributions to H, I & M */
256 store_H = (double *)space(sizeof(double) * (o+2));
257 /* unpaired between ]l,o[ */
258 store_Io = (double *)space(sizeof(double) * (o+2));
259 /* qm[p+1,i-1]*dangles_po */
260 store_M_qm_o = (double *)space(sizeof(double) * (n+1));
262 for (p=o-TURN-1; p>=1; p--) {
263 /* construction of partition function of segment [p,o], given that
264 an unpaired region [i,i+w[ exists within [p,o] */
270 /*hairpin contribution*/
271 if (((type==3)||(type==4))&&no_closingGU)
274 temp = prpr[po] * exp_E_Hairpin(u, type, S1[p+1], S1[o-1], sequence+p-1, Pf) * scale[u+2];
275 /* all H contribs are collect for the longest unpaired region */
278 /* interior loops with interior pair k,l and an unpaired region of
279 length w between p and k || l and o*/
280 for (k=p+1; k<=MIN2(p+MAXLOOP+1,o-TURN-2); k++) {
283 for (l=MAX2(k+TURN+1,o-1-MAXLOOP+u1); l<o; l++) {
286 if((l+1) < o) store_Io[l+1] += sum_l;
290 type_2 = rtype[type_2];
291 temp = prpr[po] * qb[kl] * exp_E_IntLoop(u1, o-l-1, type, type_2, S1[p+1], S1[o-1], S1[k-1], S1[l+1], Pf) *scale[u1+o-l+1];
292 if((l+1) < o) store_Io[l+1] += temp; /* unpaired region between ]l,o[ */
294 } /* end of if pair(k,l) */
296 /* unpaired in region ]p,k[ */
297 for(i=p+1;i <= k-1;i++)
298 store_I2o[i][MIN2(w-1,k-i-1)] += sum_l;
300 } /*end of if(type) test for bp (p,o) */
302 /* multiple stem loop contribution
303 calculate qm2[my_iindx[i]-j] in the course of the calculation
304 of the multiple stem loop contribution:
305 advantage: you save memory:
306 instead of a (n+1)*n array for qqm2 you only need 2*n arrays
307 disadvantage: you have to use two times the op-loop for the full
308 multiloop contribution
309 first op-loop: index o goes from 1...n and
310 index p from o-TURN-1 ... 1
311 second op-loop: index o goes from n...1 and
312 index p from o+TURN+1 ... n !!
313 HERE index o goes from 1...n and index p o-TURN-1 ... 1 ,
314 we calculate the contributions to multiple stem loop
315 where exp(i+w-1-p)*(qqm2 values between i+w and o-1)
316 AND qm[iindex[p+1]-(i-1)]*exp(beta*w)*qm[iindex[i+w]-(o-1)]
317 you have to recalculate of qqm matrix containing final stem
318 contributions to multiple loop partition function
322 qqm[p] := (contribution with exact one loop in region (p,o)*/
323 qqm[p] = qqm1[p] * expMLbase[1];
325 qbt1 = qb[po] * exp_E_MLstem(type, (p>1) ? S1[p-1] : -1, (o<n) ? S1[o+1] : -1, Pf);
327 /* reverse dangles for prpr[po]*... */
330 temp = prpr[po] * exp_E_MLstem(tt, S1[o-1], S1[p+1], Pf) * scale[2] * Pf->expMLclosing;
331 for(i=p+1; i < o; i++) {
332 int p1i = (p+1) < (i-1) ? my_iindx[p+1]-(i-1) : 0;
333 /*unpaired region expMLbase[i-p] left of structured
335 /* @expMLbase: note distance of i-p == i-(p+1)+1 */
336 store_M_mlbase[my_iindx[p+1]-i] += expMLbase[i-p] * temp * qq_1m2[i+1];
337 /* structured region qm[p1i] left of unpaired region */
338 /* contribition for unpaired region is added after the p-loop */
339 store_M_qm_o[i] += qm[p1i] * temp;
340 } /*end of for i ... */
343 for(tqm2 = 0., i=p+1; i < o; i++)
344 tqm2 += qm[my_iindx[p]-i] * qqm[i+1];
346 /* qqm2[p] contrib with at least 2 loops in region (p,o) */
348 } /* end for (p=..) */
350 for(sum_h = 0., i=1; i < o; i++) {
353 max_v = MIN2(w-1,o-i-1);
354 for(v=max_v; v >= 0; v--){
356 pu_test->H[i][v] += sum_h;/* store_H[i][v] + store_H[i][max_v]; */
357 /* Interior loops: unpaired region between ]l,o[ calculated here !*/
358 /* unpaired region between ]p,k[ collected after after o-loop */
359 if(v <= MIN2(max_v,MAXLOOP)) {
360 pu_test->I[i][v] += store_Io[i]; /* ]l,o[ */
363 /* unpaired region [i,v] between structured regions ]p,i[ and ]v,o[. */
364 /* store_M_qm_o[i] = part. funct over all structured regions ]p,i[ */
365 vo = (i+v+1) <= (o-1) ? my_iindx[i+v+1]-(o-1): 0;
366 pu_test->M[i][v] += store_M_qm_o[i]*expMLbase[v+1]*qm[vo];
369 tmp = qqm1; qqm1=qqm; qqm=tmp;
370 tmp = qqm2; qqm2=qq_1m2; qq_1m2=tmp;
375 }/* end for (o=..) */
377 for(i=1; i < n; i++) {
381 max_v = MIN2(w-1,n-i);
382 for(v=n; v >=0; v--) {
383 if(v <= MIN2(max_v,MAXLOOP)) {
384 /* all unpaired regions [i,v] between p and k in interior loops */
385 /* notice v runs from max_v -> 0, sum_iv sums all int. l. contribs */
386 /* for each x, v < x =< max_v, since they contribute to [i,v] */
387 sum_iv += store_I2o[i][v];
388 pu_test->I[i][v] += sum_iv;
390 /* all unpaired region [i,v] for a fixed v, given that */
391 /* region ]v,o[ contains at least 2 structures qq_1m2[v+1]; */
393 sum_M[v] += store_M_mlbase[my_iindx[i]-v];
395 pu_test->M[i][v-i] += sum_M[v];
406 for (i=1; i<=n; i++) {
407 /* set auxillary arrays to 0 */
408 qqm[i] = qqm1[i] = qqm2[i] = qq_1m2[i] = 0;
411 /* 2. exterior bp (p,o) encloses unpaired region [i,j]
412 HERE index o goes from n...1 and index p from o+TURN+1 ... n,
413 that is, we add the one multiloop contribution that we
414 could not calculate before */
416 /* is free'ing plus allocating faster than looping over all entries an setting them to 0? */
418 free(store_M_mlbase);
419 store_M_mlbase = (double *) space(sizeof(double) * (size + 1));
421 /* this should be the fastest way to set everything to 0 */
422 memset(store_M_mlbase, 0, sizeof(double) * (size + 1));
425 for (o=n-TURN-1;o>=1; o--) {
426 for (p=o+TURN+1; p<=n; p++) {
429 /* recalculate of qqm matrix containing final stem
430 contributions to multiple loop partition function
431 from segment [o,p] */
432 qqm[p] = qqm1[p] * expMLbase[1];
435 qbt1 *= exp_E_MLstem(type, (o>1) ? S1[o-1] : -1, (p<n) ? S1[p+1] : -1, Pf);
437 /* revers dangles for prpr[po]... */
440 temp = prpr[po]*exp_E_MLstem(tt, S1[p-1], S1[o+1], Pf) * Pf->expMLclosing * scale[2];
443 for(i=o+1; i < p; i++) {
445 tqm2+=qqm[i]*qm[my_iindx[i+1]-p];
450 /* structured region qq_1m2[i-1] left of unpaired r. expMLbase[p-i]*/
451 /* @expMLbase: note distance of p-i == p+1-i+1 */
452 store_M_mlbase[my_iindx[i]-p+1] += qq_1m2[i-1]*expMLbase[p-i]*temp;
454 }/*end of for i ....*/
456 }/* end for (p=..) */
457 tmp = qqm1; qqm1=qqm; qqm=tmp;
458 tmp = qqm2; qqm2=qq_1m2; qq_1m2=tmp;
459 }/* end for (o=..) */
460 /* now collect the missing multiloop contributions */
461 for(i=0;i<=n;i++) { sum_M[i]=0.; }
463 int v_max = MIN2(w-1,n-i);
465 sum_M[i] += store_M_mlbase[my_iindx[i]-v];
466 if ((v-i <= v_max) ) {
467 pu_test->M[i][v-i] += sum_M[i];
472 /* 1. region [i,j] exterior to all loops */
474 for (i=1; i<=n; i++) {
476 for(j=i; j<MIN2(i+w,n+1);j++){
478 temp=q1k[i-1]*1*scale[j-i+1]*qln[j+1]/q1k[n];
479 pu_test->E[i][j-i]+=temp;
485 free(store_M_mlbase);
491 PRIVATE void get_interact_arrays(unsigned int n1,
503 *p_c_S = (double **)space(sizeof(double *)*(n1+1));
505 for (i=1; i<=n1; i++){
506 pc_size = MIN2((w + incr5 + incr3), (int)n1);
507 (*p_c_S)[i] = (double *)space(sizeof(double) * (pc_size + 1));
508 for (j=0; j < pc_size; j++)
509 (*p_c_S)[i][j] = p_c->H[i][j] + p_c->I[i][j] + p_c->M[i][j] + p_c->E[i][j];
513 (*p_c2_S) = (double **)space(sizeof(double *) * (n2 + 1));
514 for (i=1; i<=n2; i++){
515 pc_size = MIN2(w, (int)n2);
516 (*p_c2_S)[i] = (double *)space(sizeof(double) * (pc_size + 2));
517 for (j=0; j < pc_size; j++)
518 (*p_c2_S)[i][j] = p_c2->H[i][j] + p_c2->I[i][j] + p_c2->M[i][j] + p_c2->E[i][j];
523 /*------------------------------------------------------------------------*/
524 /* s1 is the longer seq */
525 PUBLIC interact *pf_interact( const char *s1,
534 int i, j, k,l,n1,n2,add_i5,add_i3,i_max,k_max, pc_size;
535 double temp, Z, rev_d, E, Z2,**p_c_S, **p_c2_S, int_scale;
536 FLT_OR_DBL ****qint_4, **qint_ik;
537 /* PRIVATE double **pint; array for pf_up() output */
539 double G_min, G_is,Gi_min;
540 int gi,gj,gk,gl,ci,cj,ck,cl,prev_k,prev_l;
542 double Z_int, temp_int, temppfs;
543 double const_scale, const_T;
544 constrain *cc = NULL; /* constrains for cofolding */
545 char *Seq, *i_long,*i_short,*pos=NULL; /* short seq appended to long one */
546 /* int ***pu_jl; */ /* positions of interaction in the short RNA */
548 G_min = G_is = Gi_min = 100.0;
549 gi = gj = gk = gl = ci = cj = ck = cl = 0;
551 n1 = (int) strlen(s1);
552 n2 = (int) strlen(s2);
556 i_long = (char *) space (sizeof(char)*(n1+1));
557 i_short = (char *) space (sizeof(char)*(n2+1));
558 Seq = (char *) space (sizeof(char)*(n1+n2+2));
563 set_encoded_seq(s1, &S, &S1);
564 set_encoded_seq(s2, &SS, &SS2);
566 cc = get_ptypes(Seq,cstruc);
568 get_interact_arrays(n1, n2, p_c, p_c2, w, incr5, incr3, &p_c_S, &p_c2_S);
570 /*array for pf_up() output */
571 Int = (interact *) space(sizeof(interact)*1);
572 Int->Pi = (double *) space(sizeof(double)*(n1+2));
573 Int->Gi = (double *) space(sizeof(double)*(n1+2));
575 /* use a different scaling for pf_interact*/
576 scale_int(s2, s1, &int_scale);
578 /* set the global scale array and the global variable pf_scale to the
579 values used to scale the interaction, keep their former values !! */
581 pf_scale = int_scale;
583 /* in order to scale expLoopEnergy correctly call*/
584 scale_stru_pf_params((unsigned) n1);
586 qint_ik = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *) * (n1+1));
587 for (i=1; i<=n1; i++) {
588 qint_ik[i] = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL) * (n1+1));
591 int_ik = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *) * (n1+1));
592 for (i=1; i<=n1; i++) {
593 int_ik[i] = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL) * (n1+1));
596 /* Gint = ( -log(int_ik[gk][gi])-( ((int) w/2)*log(pf_scale)) )*((Pf->temperature+K0)*GASCONST/1000.0); */
597 const_scale = ((int) w/2)*log(pf_scale);
598 const_T = (Pf->kT/1000.0);
600 /* static short *S~S1, *S1~SS1, *SS~S2, *SS2; */
601 for (i=0; i<=n1; i++) {
602 Int->Pi[i]=Int->Gi[i]=0.;
607 if ( fold_constrained && cstruc != NULL) {
608 pos = strchr(cstruc,'|');
611 /* long seq & short seq
612 .........||..|||||....&....||||... w = maximal interaction length
614 strncpy(i_long,cstruc,n1);
616 strncpy(i_short,&cstruc[n1],n2);
618 pos = strchr(i_long,'|');
619 if(pos) ck = (int) (pos-i_long)+1; /* k */
620 pos = strrchr(i_long,'|');
621 if(pos) ci = (int) (pos-i_long)+1; /* i */
622 pos = strrchr(i_short,'|');
623 if(pos) cl = (int) (pos-i_short)+1; /* l */
624 pos = strchr(i_short,'|');
625 if(pos) cj = (int) (pos-i_short)+1; /* j */
627 if(ck > 0 && ci > 0 && ci-ck+1 > w) {
628 fprintf(stderr, "distance between constrains in longer seq, %d, larger than -w = %d",ci-ck+1,w);
629 nrerror("pf_interact: could not satisfy all constraints");
631 if(cj > 0 && cl > 0 && cl-cj+1 > w) {
632 fprintf(stderr, "distance between constrains in shorter seq, %d, larger than -w = %d",cl-cj+1,w);
633 nrerror("pf_interact: could not satisfy all constraints");
637 } else if ( fold_constrained && cstruc == NULL) {
638 nrerror("option -C selected, but no constrained structure given\n");
640 if(fold_constrained) pos = strchr(cstruc,'|');
642 /* qint_4[i][j][k][l] contribution that region (k-i) in seq1 (l=n1)
643 is paired to region (l-j) in seq 2(l=n2) that is
644 a region closed by bp k-l and bp i-j */
645 qint_4 = (FLT_OR_DBL ****) space(sizeof(FLT_OR_DBL ***) * (n1+1));
647 /* qint_4[i][j][k][l] */
648 for (i=1; i<=n1; i++) {
651 if(fold_constrained && pos && ci) end_k= MAX2(i-w, ci-w);
652 /* '|' constrains for long sequence: index i from 1 to n1 (5' to 3')*/
653 /* interaction has to include 3' most '|' constrain, ci */
654 if(fold_constrained && pos && ci && i==1 && i<ci)
655 i= ci-w+1 > 1 ? ci-w+1 : 1;
656 /* interaction has to include 5' most '|' constrain, ck*/
657 if(fold_constrained && pos && ck && i > ck+w-1) break;
659 /* note: qint_4[i] will be freed before we allocate qint_4[i+1] */
660 qint_4[i] = (FLT_OR_DBL ***) space(sizeof(FLT_OR_DBL **) * (n2+1));
661 for (j=n2; j>0; j--) {
662 qint_4[i][j] = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL*) * (w+1));
663 for (k=0; k<=w; k++) {
664 qint_4[i][j][k] = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL) * (w+1));
669 for (j=n2; j>0; j--) {
670 int type, type2,end_l;
672 if(fold_constrained && pos && ci) end_l= MIN2(cj+w,j+w);
673 /* '|' constrains for short sequence: index j from n2 to 1 (3' to 5')*/
674 /* interaction has to include 5' most '|' constrain, cj */
675 if(fold_constrained && pos && cj && j==n2 && j>cj)
676 j = cj+w-1 > n2 ? n2 : cj+w-1;
677 /* interaction has to include 3' most '|' constrain, cl*/
678 if(fold_constrained && pos && cl && j < cl-w+1) break;
679 type = cc->ptype[cc->indx[i]-(n1+j)];
680 qint_4[i][j][0][0] = type ? Pf->expDuplexInit : 0;
683 qint_4[i][j][0][0] *= exp_E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, Pf);
685 rev_d = exp_E_ExtLoop(rtype[type], (j>1) ? SS2[j-1] : -1, (i<n1) ? S1[i+1] : -1, Pf);
687 /* add inc5 and incr3 */
688 if((i-incr5) > 0 ) add_i5=i-incr5;
691 pc_size = MIN2((w+incr3+incr5),n1);
692 if(incr3 < pc_size) add_i3=incr3;
693 else add_i3=pc_size-1;
695 /* only one bp (no interior loop) */
696 if(p_c2 == NULL) {/* consider only structure of longer seq. */
697 qint_ik[i][i]+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*scale[((int) w/2)];
698 Z+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*scale[((int) w/2)];
699 } else {/* consider structures of both seqs. */
700 qint_ik[i][i]+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*p_c2_S[j][0]*scale[((int) w/2)];
701 Z+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*p_c2_S[j][0]*scale[((int) w/2)];
705 /* check deltaG_ges = deltaG_int + deltaG_unstr; */
706 int_ik[i][i]+=qint_4[i][j][0][0]*rev_d*scale[((int) w/2)];
707 Z_int+=qint_4[i][j][0][0]*rev_d*scale[((int) w/2)];
712 for (k=i-1; k>end_k && k>0; k--) {
713 if (fold_constrained && pos && cstruc[k-1] == '|' && k > prev_k)
715 for (l=j+1; l< end_l && l<=n2; l++) {
717 double scalew, tt, intt;
719 type2 = cc->ptype[cc->indx[k]-(n1+l)];
720 /* '|' : l HAS TO be paired: not pair (k,x) where x>l allowed */
721 if(fold_constrained && pos && cstruc[n1+l-1] == '|' && l < prev_l)
723 if(fold_constrained && pos && (k<=ck || i>=ci) && !type2) continue;
724 if(fold_constrained && pos && ((cstruc[k-1] == '|') || (cstruc[n1+l-1] == '|')) && !type2) break;
726 if (!type2) continue;
727 /* to save memory keep only qint_4[i-w...i][][][] in memory
728 use indices qint_4[i][j][a={0,1,...,w-1}][b={0,1,...,w-1}] */
729 a=i-k;/* k -> a from 1...w-1*/
730 b=l-j;/* l -> b from 1...w-1 */
732 /* scale everything to w/2 */
735 scalew = ( scale[isw - (a+b)] );
736 } else if ( (a+b) > isw ) {
737 scalew = 1/( scale[(a+b) - isw] );
742 if (i-k+l-j-2<=MAXLOOP) {
743 if(k >= prev_k && l <= prev_l) { /* don't violate constrains */
744 E = exp_E_IntLoop(i-k-1,l-j-1, type2, rtype[type],
745 S1[k+1], SS2[l-1], S1[i-1], SS2[j+1], Pf) *
746 scale[i-k+l-j]; /* add *scale[u1+u2+2] */
748 qint_4[i][j][a][b] += ( qint_4[k][l][0][0]*E);
750 /* use ia and ib to go from a....w-1 and from b....w-1 */
752 while((a+ia)<w && i-(a+ia)>=1 && (b+ib)<w && (j+b+ib)<=n2) {
755 qint_4[i][j][a+ia][b+ib] += qint_4[k][l][ia][ib]*E;
758 while(a+iaa<w && i-(a+iaa)>=1) {
759 qint_4[i][j][a+iaa][b+ib] += qint_4[k][l][iaa][ib]*E;
764 while( (b+ibb)<w && (j+b+ibb)<=n2 ) {
765 qint_4[i][j][a+ia][b+ibb] += qint_4[k][l][ia][ibb]*E;
773 /* '|' constrain in long sequence */
774 /* collect interactions starting before 5' most '|' constrain */
775 if ( fold_constrained && pos && ci && i < ci) continue;
776 /* collect interactions ending after 3' most '|' constrain*/
777 if ( fold_constrained && pos && ck && k > ck) continue;
778 /* '|' constrain in short sequence */
779 /* collect interactions starting before 5' most '|' constrain */
780 if ( fold_constrained && pos && cj && j > cj) continue;
781 /* collect interactions ending after 3' most '|' constrain*/
782 if ( fold_constrained && pos && cl && l < cl) continue;
784 /* scale everything to w/2*/
785 /* qint_ik[k][i] all interactions where k and i both are paired */
786 /* substract incr5 from k */
787 if(k-incr5 > 0) add_i5=k-incr5;
790 pc_size = MIN2((w+incr3+incr5),n1);
791 if(i-k+incr3 < pc_size) add_i3=i-k+incr3;
792 else add_i3=pc_size-1;
794 if(p_c2 == NULL) {/* consider only structure of longer seq. */
795 tt = qint_4[i][j][a][b]*p_c_S[add_i5][add_i3]*scalew*rev_d;
796 } else { /* consider structures of both seqs. */
797 tt = qint_4[i][j][a][b]*p_c_S[add_i5][add_i3]*p_c2_S[j][b]*scalew*rev_d;
802 /* check deltaG_ges = deltaG_int + deltaG_unstr; */
803 intt = qint_4[i][j][a][b]*scalew*rev_d;
806 G_is = (-log(tt)-const_scale)*(const_T);
807 if (G_is < G_min || EQUAL(G_is,G_min)) {
809 Gi_min =(-log(intt)-const_scale)*(const_T);
822 /* free qint_4 values not needed any more */
826 if (fold_constrained && pos && ci && i-w < ci-w+1) continue;
827 if (fold_constrained && pos && ci) bla = MAX2(ci-w+1,i-w);
828 for (j=n2; j>0; j--) {
829 for (k=0; k<=w; k++){
830 free(qint_4[bla][j][k]);
832 free(qint_4[bla][j]);
842 for (i=1; i<=n1; i++) {
843 for (k=i; k<=n1 && k<i+w; k++) {
846 /* Int->Pi[l]: prob that position l is within a paired region */
847 /* qint_ik[i][k] as well as Z are scaled to scale[((int) w/2) */
848 Int->Pi[l]+=qint_ik[i][k]/Z;
849 /* Int->Gi[l]: minimal delta G at position [l] */
850 Int->Gi[l]=MIN2(Int->Gi[l],
851 ( -log(qint_ik[i][k])-( ((int) w/2)*log(pf_scale)) )*
860 if (fold_constrained && pos && ci) {
861 /* a break in the k loop might result in unfreed values */
862 start_i = ci-w+1 < n1-w+1 ? ci-w+1 : n1-w+1;
863 start_i = start_i > 0 ? start_i : 1;
865 end_i = ck+w-1 > n1 ? n1 : ck+w-1;
867 for (i=start_i; i<=end_i; i++) {
868 if(qint_4[i] == NULL ) continue;
869 for (j=n2; j>0; j--) {
870 for (k=0; k<=w; k++) {
871 free(qint_4[i][j][k]);
882 if (fold_constrained && pos) {
883 start_i = ci-w+1 > 0 ? ci-w+1 : 1;
884 end_i = ck+w-1 > n1 ? n1 : ck+w-1;
887 for (i=start_i; i<=end_i; i++) {
888 for (j=n2; j>0; j--) {
889 for (k=0; k<=w; k++) {
890 free(qint_4[i][j][k]);
898 if(fold_constrained && (gi==0 || gk==0 || gl==0 || gj==0)) {
899 nrerror("pf_interact: could not satisfy all constraints");
901 /* fill structure interact */
908 Int->Gikjl_wo = Gi_min;
913 for (i=1; i<=n1; i++) {
917 for (i=1; i<=n1; i++) {
922 /* reset the global variables pf_scale and scale to their original values */
923 pf_scale = temppfs;/* reset pf_scale */
924 scale_stru_pf_params((unsigned) n1);/* reset the scale array */
925 free_pf_arrays(); /* for arrays for pf_fold(...) */
927 if(expMLbase != NULL) {
935 for (i=1; i<=n1; i++) {
940 for (i=1; i<=n2; i++) {
951 /*------------------------------------------------------------------------*/
952 /* use an extra scale for pf_interact, here sl is the longer sequence */
953 PRIVATE void scale_int(const char *s, const char *sl, double *sc_int){
961 expMLbase = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(nl+1));
962 scale = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*((nl+1)*2));
964 /* use RNA duplex to get a realistic estimate for the best possible
965 interaction energy between the short RNA s and its target sl */
966 mfe = duplexfold(s,sl);
968 kT = Pf->kT/1000.0; /* in Kcal */
970 /* sc_int is similar to pf_scale: i.e. one time the scale */
971 *sc_int = exp(-(mfe.energy)/kT/n);
973 /* free the structure returned by duplexfold */
977 /*----------------------------------------------------------------------*/
978 /* init_pf_two(n) :gets the arrays, that you need, from part_func.c */
979 /* get_pf_arrays(&S, &S1, &ptype, &qb, &qm, &q1k, &qln);*/
980 /* init_pf_fold(), update_pf_params, encode_char(), make_ptypes() are called by pf_fold() */
981 PRIVATE void init_pf_two(int length){
983 nonstandard_arithmetic();
991 /* gets the arrays, that we need, from part_func.c */
992 if(!get_pf_arrays(&S, &S1, &ptype, &qb, &qm, &q1k, &qln))
993 nrerror("init_pf_two: pf_fold() has to be called before calling pf_unstru()\n");
994 /* get a pointer to the base pair probs */
995 probs = export_bppm();
997 scale_stru_pf_params((unsigned) length);
1000 if(init_temp != Pf->temperature)
1001 nrerror("init_pf_two: inconsistency with temperature");
1004 PRIVATE void get_up_arrays(unsigned int length){
1005 unsigned int l1 = length + 1;
1006 unsigned int l2 = length + 2;
1007 prpr = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * ((l1*l2)>>1));
1008 expMLbase = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * l1);
1009 scale = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * l1);
1010 qqm2 = (double *) space(sizeof(double) * l2);
1011 qq_1m2 = (double *) space(sizeof(double) * l2);
1012 qqm = (double *) space(sizeof(double) * l2);
1013 qqm1 = (double *) space(sizeof(double) * l2);
1014 my_iindx = get_iindx(length);
1017 PRIVATE void free_up_arrays(void){
1018 if(prpr != NULL){ free(prpr); prpr = NULL;}
1019 if(expMLbase != NULL){ free(expMLbase); expMLbase = NULL;}
1020 if(scale != NULL){ free(scale); scale = NULL;}
1021 if(qqm != NULL){ free(qqm); qqm = NULL;}
1022 if(qqm1 != NULL){ free(qqm1); qqm1 = NULL;}
1023 if(qqm2 != NULL){ free(qqm2); qqm2 = NULL;}
1024 if(qq_1m2 != NULL){ free(qq_1m2); qq_1m2 = NULL;}
1025 if(my_iindx != NULL){ free(my_iindx); my_iindx = NULL;}
1028 PUBLIC void free_interact(interact *pin) {
1029 if(S != NULL && pin != NULL){
1033 if(S1 != NULL && pin != NULL){
1044 /*---------------------------------------------------------------------------*/
1046 PRIVATE void encode_seq(const char *s1, const char *s2) {
1050 /* S and S1 are freed by free_pf_arrays(); ! */
1051 S = (short *) space(sizeof(short)*(l+1));
1052 S1= (short *) space(sizeof(short)*(l+1));
1053 /* S1 exists only for the special X K and I bases and energy_set!=0 */
1055 for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
1056 S[i]= (short) encode_char(toupper(s1[i-1]));
1057 S1[i] = alias[S[i]]; /* for mismatches of nostandard bases */
1061 /* SS2 exists only for the special X K and I bases and energy_set!=0 */
1063 for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
1064 SS[i]= (short) encode_char(toupper(s2[i-1]));
1065 SS2[i] = alias[SS[i]]; /* for mismatches of nostandard bases */
1070 /*-------------------------------------------------------------------------*/
1071 /* scale energy parameters and pre-calculate Boltzmann weights:
1072 most of this is done in structure Pf see params.c,h (function:
1073 get_scaled_pf_parameters(), only arrays scale and expMLbase are handled here*/
1074 PRIVATE void scale_stru_pf_params(unsigned int length)
1080 /* Do this only at the first call for get_scaled_pf_parameters()
1081 and/or if temperature has changed*/
1082 if(init_temp != temperature) {
1084 Pf=get_scaled_pf_parameters();
1087 init_temp = Pf->temperature;
1089 kT = Pf->kT; /* kT in cal/mol */
1091 /* scaling factors (to avoid overflows) */
1092 if (pf_scale == -1) { /* mean energy for random sequences: 184.3*length cal */
1093 pf_scale = exp(-(-185+(Pf->temperature-37.)*7.27)/kT);
1094 if (pf_scale<1) pf_scale=1;
1097 scale[1] = 1./pf_scale;
1099 expMLbase[1] = Pf->expMLbase/pf_scale;
1100 for (i=2; i<=length; i++) {
1101 scale[i] = scale[i/2]*scale[i-(i/2)];
1102 expMLbase[i] = pow(Pf->expMLbase, (double)i) * scale[i];
1105 /*-------------------------------------------------------------------------*/
1106 /* make a results structure containing all u-values & the header */
1107 PUBLIC pu_out *get_u_vals(pu_contrib *p_c, int **unpaired_values, char *select_contrib) {
1108 int i, j, k, l, g,ws,num_u_vals,unstr,count,contribs,size,w,len;
1110 int off_S, off_E, off_H, off_I, off_M;
1111 /* double **p_cont,**p_cont_sh, dG_u; p_u AND its contributions */
1116 /* number of different -u values */
1117 for (num_u_vals = 0, i = 1; i <= unpaired_values[0][0]; i++) {
1118 j = unpaired_values[i][0];
1119 do num_u_vals++; while(++j <= unpaired_values[i][1]);
1121 /* check which contributions ([-c "SHIME"] ) are desired by the user,
1122 set the offset for each contribution */
1124 S = E = H = I = M = 0;
1125 off_S = off_E = off_H = off_I = off_M = 0;
1126 if(strchr(select_contrib, 'S')) {
1131 if(strchr(select_contrib, 'E')) {
1136 if(strchr(select_contrib, 'H')) {
1141 if(strchr(select_contrib, 'I')) {
1146 if(strchr(select_contrib, 'M')) {
1153 nrerror("get_u_vals: error with contribs!");
1155 /* allocate the results structure */
1156 u_results = (pu_out *) space(1*sizeof(pu_out));
1157 u_results->len = len; /* sequence length */
1158 /*num_u_vals differnet -u values, contribs [-c "SHIME"] */
1159 u_results->u_vals = num_u_vals;
1160 u_results->contribs = contribs;
1161 /* add 1 column for position within the sequence and
1162 add 1 column for the free energy of interaction values */
1163 /* header e.g. u3I (contribution for u3 interior loops */
1164 size = 1 + (num_u_vals*contribs) + 1;
1165 u_results->header = (char **) space((size+1)*sizeof(char*));
1166 for(i=0;i<(size+1);i++){
1167 u_results->header[i] = (char *) space(10*sizeof(char));
1169 /* different free energies for all -u and -c combinations */
1170 u_results->u_values = (double**) space((size+1) *sizeof(double*));
1171 for(i=0;i<(size+1);i++){
1172 /* position within the sequence */
1173 u_results->u_values[i] = (double*) space((len+3)*sizeof(double));
1175 /* write the position within the sequence in the u_results array
1177 sprintf(u_results->header[0],"pos");
1178 for(i=0;i<=len;i++){
1179 /* add the position*/
1180 u_results->u_values[0][i] = i;
1182 /* go over the different -u values, u_vals[] listy of different -u values*/
1183 for (count = k = 1; k <= unpaired_values[0][0]; k++) {
1184 l = unpaired_values[k][0];
1186 int offset; /* offset for the respective -u value (depents on the number
1187 of the -u value and on the numbers of contribs */
1189 offset = ((count - 1) * contribs) + 1; /* first colum is the position */
1190 /* set the current value of -u : here we call it w */
1191 w = l; /* set w to the actual -u value */
1192 if(w > len) break; /* corr caro */
1193 /* make the header - look which contribitions are wanted */
1194 if(S) sprintf(u_results->header[offset+off_S],"u%dS",w);
1195 if(E) sprintf(u_results->header[offset+off_E],"u%dE",w);
1196 if(H) sprintf(u_results->header[offset+off_H],"u%dH",w);
1197 if(I) sprintf(u_results->header[offset+off_I],"u%dI",w);
1198 if(M) sprintf(u_results->header[offset+off_M],"u%dM",w);
1201 for (i=1; i<=len; i++) { /* for each position */
1202 /* w goes form j to i (intervall end at i) */
1203 for (j=i; j < MIN2((i+w),len+1); j++) { /* for each -u value < w
1204 this is not necessay ->
1205 calculate j from i and w
1208 /* if (j-i+1) == w we have the -u = w value wanted */
1209 if( (j-i+1) == w && i+w-1 <= len) {
1210 blubb = p_c->H[i][j-i]+p_c->I[i][j-i]+p_c->M[i][j-i]+p_c->E[i][j-i];
1212 /* printf("len %d blubb %.3f \n",len, blubb); */
1213 if(S) u_results->u_values[offset+off_S][i+w-1]+=blubb;
1214 if(E) u_results->u_values[offset+off_E][i+w-1]+=p_c->E[i][j-i];
1215 if(H) u_results->u_values[offset+off_H][i+w-1]+=p_c->H[i][j-i];
1216 if(I) u_results->u_values[offset+off_I][i+w-1]+=p_c->I[i][j-i];
1217 if(M) u_results->u_values[offset+off_M][i+w-1]+=p_c->M[i][j-i];
1220 if(i<w && (j-i+1) != w && i+w-1 > len && i+w-1 < len+3) {
1221 if(S) u_results->u_values[offset+off_S][i+w-1]=-1;
1222 if(E) u_results->u_values[offset+off_E][i+w-1]=-1;
1223 if(H) u_results->u_values[offset+off_H][i+w-1]=-1;
1224 if(I) u_results->u_values[offset+off_I][i+w-1]=-1;
1225 if(M) u_results->u_values[offset+off_M][i+w-1]=-1;
1229 } else return(NULL); /* error */
1231 } while(++l <= unpaired_values[k][1]);
1233 return(u_results); /*success*/
1235 /* plot the results structure */
1236 /* when plotting the results for the target seq we add a header */
1237 /* when plotting the results for the interaction partner u want no header,
1238 set s1 to NULL to avoid plotting the header */
1239 /* currently we plot the free energies to a file: the probability of
1240 being unpaired for region [i,j], p_u[i,j], is related to the free
1241 energy to open region [i,j], dG_u[i,j] by:
1242 dG_u[i,j] = -log(p_u[i,j])*(temperature+K0)*GASCONST/1000.0; */
1243 PUBLIC int plot_free_pu_out(pu_out* res, interact *pint, char *ofile, char *head) {
1246 char nan[4], *time, dg[11];
1249 wastl = fopen(ofile,"a");
1251 fprintf(stderr, "p_cont: can't open %s for Up_plot\n", ofile);
1256 /* printf("T=%.16f \n(temperature+K0)*GASCONST/1000.0 = %.16f\n",temperature,(temperature+K0)*GASCONST/1000.0); */
1258 /* write the header of the output file: */
1259 /* # timestamp commandlineaufruf */
1260 /* # length and name of first sequence (target)
1262 /* # length and name of second sequence (interaction partner) */
1264 /* the next line is the output for the target: colums
1265 position in target | dG_unpaired values for target | interaction energy */
1266 /* # pos u1S u1H dg */
1267 /* values for target */
1268 /* if -b was choosen: the next lines are the dG_unpaired values for
1269 the interaction partner */
1271 /* values for the interaction partner */
1273 /* print header, if nh is zerro */
1275 time = time_stamp();
1276 fprintf(wastl,"# %s\n", time);
1277 fprintf(wastl,"%s\n",head);
1279 fprintf(wastl,"# ");
1280 /* } else { fprintf(wastl," "); } close if before */
1282 size = res->u_vals * res->contribs;
1287 for(i=0;i<=len; i++) {
1288 for(s=0;s<=size+1;s++) { /* that is for different contribution */
1289 if ( i== 0 && s > size && pint != NULL)
1290 fprintf(wastl,"%8s ",dg);
1292 if(s>0 && s<=size) {
1293 if(res->u_values[s][i] > 0.0) {
1294 dG_u = -log(res->u_values[s][i])*kT/1000.0;
1295 fprintf(wastl,"%8.3f ",dG_u);
1296 } else { /* no p_u value was defined print nan*/
1297 fprintf(wastl,"%8s ",nan);
1300 } else if (s > size && pint != NULL) {
1301 fprintf(wastl,"%8.3f ",pint->Gi[i]);
1302 } else if (s == 0) {
1303 fprintf(wastl,"%8.0f ",res->u_values[s][i]);
1307 fprintf(wastl,"%8s ",res->header[s]);
1309 fprintf(wastl,"%7s ",res->header[s]);
1313 fprintf(wastl,"\n");
1316 /*free pu_out* res */
1318 for(i=0;i<=(size+2);i++) {
1319 free(res->u_values[i]);
1320 free(res->header[i]);
1322 free(res->u_values);
1328 return(1); /* success */
1331 PUBLIC int Up_plot(pu_contrib *p_c, pu_contrib *p_c_sh, interact *pint, char *ofile, int **unpaired_values, char *select_contrib, char *head, unsigned int mode) {
1334 /* check what case we have */
1336 /* upmode = 1 only one seq */
1337 /* if(p_c != NULL && pint == NULL) { */
1338 if(mode & RNA_UP_MODE_1){
1339 dada = get_u_vals(p_c,unpaired_values,select_contrib);
1340 ret = plot_free_pu_out(dada,NULL,ofile,head);
1342 /* upmode > 1 cofolding */
1343 /* } else if (p_c != NULL && pint != NULL) { */
1344 } else if(mode & RNA_UP_MODE_2) {
1345 dada = get_u_vals(p_c,unpaired_values,select_contrib);
1346 ret = plot_free_pu_out(dada,pint,ofile,head);
1348 /* upmode = 3 cofolding*/
1349 /* } else if (p_c == NULL && p_c_sh != NULL) { */
1351 if(mode & RNA_UP_MODE_3) {
1352 dada = get_u_vals(p_c,unpaired_values, select_contrib);
1353 ret = plot_free_pu_out(dada, pint, ofile, head);
1355 dada = get_u_vals(p_c_sh, unpaired_values, select_contrib);
1356 ret = plot_free_pu_out(dada,NULL,ofile, NULL);
1361 /*-------------------------------------------------------------------------*/
1362 /* copy from part_func_co.c */
1363 PRIVATE constrain *get_ptypes(char *Seq, const char *structure) {
1364 int n,i,j,k,l, length;
1368 length = strlen(Seq);
1370 con = (constrain *) space(sizeof(constrain));
1371 con->indx = (int *) space(sizeof(int)*(length+1));
1372 for (i=1; i<=length; i++) {
1373 con->indx[i] = ((length+1-i)*(length-i))/2 +length+1;
1375 con->ptype = (char *) space(sizeof(char)*((length+1)*(length+2)/2));
1377 set_encoded_seq((const char *)Seq, &s, &s1);
1380 for (k=1; k<=n-CO_TURN-1; k++)
1381 for (l=1; l<=2; l++) {
1382 int type,ntype=0,otype=0;
1383 i=k; j = i+CO_TURN+l; if (j>n) continue;
1384 type = pair[s[i]][s[j]];
1385 while ((i>=1)&&(j<=n)) {
1386 if ((i>1)&&(j<n)) ntype = pair[s[i-1]][s[j+1]];
1387 if (noLonelyPairs && (!otype) && (!ntype))
1388 type = 0; /* i.j can only form isolated pairs */
1389 con->ptype[con->indx[i]-j] = (char) type;
1396 if (fold_constrained&&(structure!=NULL)) {
1399 stack = (int *) space(sizeof(int)*(n+1));
1400 for(hx=0, j=1; j<=n; j++) {
1401 switch (structure[j-1]) {
1402 case 'x': /* can't pair */
1403 for (l=1; l<j-CO_TURN; l++) con->ptype[con->indx[l]-j] = 0;
1404 for (l=j+CO_TURN+1; l<=n; l++) con->ptype[con->indx[j]-l] = 0;
1409 case '<': /* pairs upstream */
1413 fprintf(stderr, "%s\n", structure);
1414 nrerror("1. unbalanced brackets in constraints");
1417 type = con->ptype[con->indx[i]-j];
1418 /* don't allow pairs i<k<j<l */
1419 for (k=i; k<=j; k++)
1420 for (l=j; l<=n; l++) con->ptype[con->indx[k]-l] = 0;
1421 /* don't allow pairs k<i<l<j */
1422 for (k=1; k<=i; k++)
1423 for (l=i; l<=j; l++) con->ptype[con->indx[k]-l] = 0;
1424 con->ptype[con->indx[i]-j] = (type==0)?7:type;
1425 case '>': /* pairs downstream */
1430 fprintf(stderr, "%s\n", structure);
1431 nrerror("2. unbalanced brackets in constraint string");
1439 PRIVATE void set_encoded_seq(const char *sequence, short **S, short **S1){
1441 l = strlen(sequence);
1443 *S = (short *)space(sizeof(short) * (l + 2));
1444 for(i=1; i<=l; i++) /* make numerical encoding of sequence */
1445 (*S)[i]= (short) encode_char(toupper(sequence[i-1]));
1446 (*S)[l+1] = (*S)[1];
1447 (*S)[0] = (short) l;
1449 /* S1 exists only for the special X K and I bases and energy_set!=0 */
1451 *S1 = (short *)space(sizeof(short) * (l + 2));
1452 for(i=1; i<=l; i++) /* make numerical encoding of sequence */
1453 (*S1)[i] = alias[(short) encode_char(toupper(sequence[i-1]))]; /* for mismatches of nostandard bases */
1454 /* for circular folding add first base at position n+1 and last base at position 0 in S1 */
1455 (*S1)[l+1] = (*S1)[1];
1456 (*S1)[0] = (*S1)[l];