1 /* Last changed Time-stamp: <2009-02-24 15:17:17 ivo> */
3 minimum free energy folding
4 for a set of aligned sequences
26 #include "energy_par.h"
27 #include "fold_vars.h"
32 #include "loop_energies.h"
41 #define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
42 #define NEW_NINIO 1 /* new asymetry penalty */
43 #define MAXSECTORS 500 /* dimension for a backtrack array */
44 #define LOCALITY 0. /* locality parameter for base-pairs */
46 #define MINPSCORE -2 * UNIT
49 #################################
51 #################################
53 PUBLIC double cv_fact=1.; /* should be made static to not interfere with other threads */
54 PUBLIC double nc_fact=1.; /* should be made static to not interfere with other threads */
57 #################################
59 #################################
61 PRIVATE short **S = NULL;
62 PRIVATE short **S5 = NULL; /*S5[s][i] holds next base 5' of i in sequence s*/
63 PRIVATE short **S3 = NULL; /*Sl[s][i] holds next base 3' of i in sequence s*/
64 PRIVATE char **Ss = NULL;
65 PRIVATE unsigned short **a2s = NULL;
66 PRIVATE paramT *P = NULL;
67 PRIVATE int *indx = NULL; /* index for moving in the triangle matrices c[] and fMl[]*/
68 PRIVATE int *c = NULL; /* energy array, given that i-j pair */
69 PRIVATE int *cc = NULL; /* linear array for calculating canonical structures */
70 PRIVATE int *cc1 = NULL; /* " " */
71 PRIVATE int *f5 = NULL; /* energy of 5' end */
72 PRIVATE int *fML = NULL; /* multi-loop auxiliary energy array */
73 PRIVATE int *Fmi = NULL; /* holds row i of fML (avoids jumps in memory) */
74 PRIVATE int *DMLi = NULL; /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j]) */
75 PRIVATE int *DMLi1 = NULL; /* MIN(fML[i+1,k]+fML[k+1,j]) */
76 PRIVATE int *DMLi2 = NULL; /* MIN(fML[i+2,k]+fML[k+1,j]) */
77 PRIVATE int *pscore = NULL; /* precomputed array of pair types */
78 PRIVATE int init_length = -1;
79 PRIVATE sect sector[MAXSECTORS]; /* stack of partial structures for backtracking */
80 PRIVATE bondT *base_pair2 = NULL;
81 PRIVATE int circular = 0;
82 PRIVATE int with_gquad = 0;
84 PRIVATE int *ggg = NULL; /* minimum free energies of the gquadruplexes */
85 PRIVATE char *cons_seq = NULL;
86 PRIVATE short *S_cons = NULL;
90 #pragma omp threadprivate(S, S5, S3, Ss, a2s, P, indx, c, cc, cc1, f5, fML, Fmi, DMLi, DMLi1, DMLi2,\
91 pscore, init_length, sector, base_pair2,\
92 ggg, with_gquad, cons_seq, S_cons)
97 #################################
98 # PRIVATE FUNCTION DECLARATIONS #
99 #################################
102 PRIVATE void init_alifold(int length);
103 PRIVATE void get_arrays(unsigned int size);
104 PRIVATE void make_pscores(const short *const *S, const char **AS, int n_seq, const char *structure);
105 PRIVATE int fill_arrays(const char **strings);
106 PRIVATE void backtrack(const char **strings, int s);
108 PRIVATE void energy_of_alistruct_pt(const char **sequences,short * ptable, int n_seq, int *energy);
109 PRIVATE void stack_energy_pt(int i, const char **sequences, short *ptable, int n_seq, int *energy);
110 PRIVATE int ML_Energy_pt(int i, int n_seq, short *pt);
111 PRIVATE int EL_Energy_pt(int i, int n_seq, short *pt);
113 PRIVATE void en_corr_of_loop_gquad(int i,
115 const char **sequences,
116 const char *structure,
123 #################################
124 # BEGIN OF FUNCTION DEFINITIONS #
125 #################################
128 /* unsafe function that will be replaced by a threadsafe companion in the future */
129 PRIVATE void init_alifold(int length){
132 /* Explicitly turn off dynamic threads */
136 if (length < 1) nrerror("initialize_fold: argument must be greater 0");
137 free_alifold_arrays();
138 get_arrays((unsigned) length);
139 init_length = length;
141 indx = get_indx((unsigned)length);
143 update_alifold_params();
146 PRIVATE void get_arrays(unsigned int size){
147 if(size >= (unsigned int)sqrt((double)INT_MAX))
148 nrerror("get_arrays@alifold.c: sequence length exceeds addressable range");
150 c = (int *) space(sizeof(int)*((size*(size+1))/2+2));
151 fML = (int *) space(sizeof(int)*((size*(size+1))/2+2));
152 pscore = (int *) space(sizeof(int)*((size*(size+1))/2+2));
153 f5 = (int *) space(sizeof(int)*(size+2));
154 cc = (int *) space(sizeof(int)*(size+2));
155 cc1 = (int *) space(sizeof(int)*(size+2));
156 Fmi = (int *) space(sizeof(int)*(size+1));
157 DMLi = (int *) space(sizeof(int)*(size+1));
158 DMLi1 = (int *) space(sizeof(int)*(size+1));
159 DMLi2 = (int *) space(sizeof(int)*(size+1));
160 if(base_pair2) free(base_pair2);
162 base_pair2 = (bondT *) space(sizeof(bondT)*(1+size/2+200)); /* add a guess of how many G's may be involved in a G quadruplex */
165 PUBLIC void free_alifold_arrays(void){
172 if(pscore) free(pscore);
173 if(base_pair2) free(base_pair2);
176 if(DMLi1) free(DMLi1);
177 if(DMLi2) free(DMLi2);
180 if(cons_seq) free(cons_seq);
181 if(S_cons) free(S_cons);
183 indx = c = fML = f5 = cc = cc1 = Fmi = DMLi = DMLi1 = DMLi2 = ggg = NULL;
194 PUBLIC void alloc_sequence_arrays(const char **sequences, short ***S, short ***S5, short ***S3, unsigned short ***a2s, char ***Ss, int circ){
195 unsigned int s, n_seq, length;
196 if(sequences[0] != NULL){
197 length = strlen(sequences[0]);
198 for (s=0; sequences[s] != NULL; s++);
200 *S = (short **) space((n_seq+1) * sizeof(short *));
201 *S5 = (short **) space((n_seq+1) * sizeof(short *));
202 *S3 = (short **) space((n_seq+1) * sizeof(short *));
203 *a2s = (unsigned short **) space((n_seq+1) * sizeof(unsigned short *));
204 *Ss = (char **) space((n_seq+1) * sizeof(char *));
205 for (s=0; s<n_seq; s++) {
206 if(strlen(sequences[s]) != length) nrerror("uneqal seqence lengths");
207 (*S5)[s] = (short *) space((length + 2) * sizeof(short));
208 (*S3)[s] = (short *) space((length + 2) * sizeof(short));
209 (*a2s)[s] = (unsigned short *)space((length + 2) * sizeof(unsigned short));
210 (*Ss)[s] = (char *) space((length + 2) * sizeof(char));
211 (*S)[s] = (short *) space((length + 2) * sizeof(short));
212 encode_ali_sequence(sequences[s], (*S)[s], (*S5)[s], (*S3)[s], (*Ss)[s], (*a2s)[s], circ);
216 (*a2s)[n_seq] = NULL;
220 else nrerror("alloc_sequence_arrays: no sequences in the alignment!");
223 PUBLIC void free_sequence_arrays(unsigned int n_seq, short ***S, short ***S5, short ***S3, unsigned short ***a2s, char ***Ss){
225 for (s=0; s<n_seq; s++) {
233 free(*S5); *S5 = NULL;
234 free(*S3); *S3 = NULL;
235 free(*a2s); *a2s = NULL;
236 free(*Ss); *Ss = NULL;
239 PUBLIC void update_alifold_params(void){
241 P = scale_parameters();
243 if (init_length < 0) init_length=0;
246 PUBLIC float alifold(const char **strings, char *structure){
247 int length, energy, s, n_seq;
250 length = (int) strlen(strings[0]);
253 /* always init everything since all global static variables are uninitialized when entering a thread */
254 init_alifold(length);
256 if (length>init_length) init_alifold(length);
258 if (fabs(P->temperature - temperature)>1e-6) update_alifold_params();
260 for (s=0; strings[s]!=NULL; s++);
262 with_gquad = P->model_details.gquad;
264 alloc_sequence_arrays(strings, &S, &S5, &S3, &a2s, &Ss, circular);
265 make_pscores((const short **) S, strings, n_seq, structure);
267 energy = fill_arrays((const char **)strings);
269 backtrack((const char **)strings, 0);
272 parenthesis_structure(structure, base_pair2, length);
274 letter_structure(structure, base_pair2, length);
278 * Backward compatibility:
279 * This block may be removed if deprecated functions
280 * relying on the global variable "base_pair" vanishs from within the package!
282 base_pair = base_pair2;
285 if(base_pair) free(base_pair);
286 base_pair = (bondT *)space(sizeof(bondT) * (1+length/2));
287 memcpy(base_pair, base_pair2, sizeof(bondT) * (1+length/2));
290 free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);
292 if (backtrack_type=='C')
293 return (float) c[indx[length]+1]/(n_seq*100.);
294 else if (backtrack_type=='M')
295 return (float) fML[indx[length]+1]/(n_seq*100.);
297 return (float) f5[length]/(n_seq*100.);
302 *** the actual forward recursion to fill the energy arrays
304 PRIVATE int fill_arrays(const char **strings) {
305 int i, j, k, p, q, length, energy, new_c;
306 int decomp, MLenergy, new_fML;
307 int s, n_seq, *type, type_2, tt;
309 /* count number of sequences */
310 for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
312 type = (int *) space(n_seq*sizeof(int));
313 length = strlen(strings[0]);
318 cons_seq = consensus(strings);
319 /* make g-island annotation of the consensus */
320 S_cons = encode_sequence(cons_seq, 0);
321 ggg = get_gquad_ali_matrix(S_cons, S, n_seq, P);
324 for (j=1; j<=length; j++){
325 Fmi[j]=DMLi[j]=DMLi1[j]=DMLi2[j]=INF;
326 for (i=(j>TURN?(j-TURN):1); i<j; i++) {
327 c[indx[j]+i] = fML[indx[j]+i] = INF;
331 /* begin recursions */
332 for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
333 for (j = i+TURN+1; j <= length; j++) {
334 int ij, psc, l1, maxq, minq, up, c0;
337 for (s=0; s<n_seq; s++) {
338 type[s] = pair[S[s][i]][S[s][j]];
339 if (type[s]==0) type[s]=7;
342 psc = pscore[indx[j]+i];
343 if (psc>=MINPSCORE) { /* a pair to consider */
344 int stackEnergy = INF;
345 /* hairpin ----------------------------------------------*/
348 for (new_c=s=0; s<n_seq; s++) {
349 if ((a2s[s][j-1]-a2s[s][i])<3) new_c+=600;
350 else new_c += E_Hairpin(a2s[s][j-1]-a2s[s][i],type[s],S3[s][i],S5[s][j],Ss[s]+(a2s[s][i-1]), P);
352 /*--------------------------------------------------------
353 check for elementary structures involving more than one
355 --------------------------------------------------------*/
357 for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
358 minq = j-i+p-MAXLOOP-2;
359 if (minq<p+1+TURN) minq = p+1+TURN;
360 for (q = minq; q < j; q++) {
361 if (pscore[indx[q]+p]<MINPSCORE) continue;
362 for (energy = s=0; s<n_seq; s++) {
363 type_2 = pair[S[s][q]][S[s][p]]; /* q,p not p,q! */
364 if (type_2 == 0) type_2 = 7;
365 energy += E_IntLoop(a2s[s][p-1]-a2s[s][i], a2s[s][j-1]-a2s[s][q], type[s], type_2,
367 S5[s][p], S3[s][q], P);
369 new_c = MIN2(new_c, energy + c[indx[q]+p]);
370 if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
374 /* multi-loop decomposition ------------------------*/
377 for(s=0; s<n_seq; s++){
379 decomp += E_MLstem(tt, S5[s][j], S3[s][i], P);
383 for(s=0; s<n_seq; s++){
385 decomp += E_MLstem(tt, -1, -1, P);
388 MLenergy = decomp + n_seq*P->MLclosing;
389 new_c = MIN2(new_c, MLenergy);
393 for(s=0;s<n_seq;s++){
396 decomp += P->mismatchI[tt][S3[s][i]][S5[s][j]];
398 decomp += P->TerminalAU;
400 for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
402 if(l1>MAXLOOP) break;
403 if(S_cons[p] != 3) continue;
404 minq = j - i + p - MAXLOOP - 2;
405 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
406 minq = MAX2(c0, minq);
408 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
409 maxq = MIN2(c0, maxq);
410 for(q = minq; q < maxq; q++){
411 if(S_cons[q] != 3) continue;
412 c0 = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[l1 + j - q - 1];
413 new_c = MIN2(new_c, c0);
419 if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
420 minq = j - i + p - MAXLOOP - 2;
421 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
422 minq = MAX2(c0, minq);
424 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
425 maxq = MIN2(c0, maxq);
426 for(q = minq; q < maxq; q++){
427 if(S_cons[q] != 3) continue;
428 c0 = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[j - q - 1];
429 new_c = MIN2(new_c, c0);
435 for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
437 if(l1>MAXLOOP) break;
438 if(S_cons[p] != 3) continue;
439 c0 = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[l1];
440 new_c = MIN2(new_c, c0);
444 new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
446 cc[j] = new_c - psc; /* add covariance bonnus/penalty */
448 c[ij] = cc1[j-1]+stackEnergy-psc;
452 } /* end >> if (pair) << */
455 /* done with c[i,j], now compute fML[i,j] */
456 /* free ends ? -----------------------------------------*/
458 new_fML = fML[ij+1]+n_seq*P->MLbase;
459 new_fML = MIN2(fML[indx[j-1]+i]+n_seq*P->MLbase, new_fML);
462 for (s=0; s<n_seq; s++) {
463 energy += E_MLstem(type[s], S5[s][i], S3[s][j], P);
467 for (s=0; s<n_seq; s++) {
468 energy += E_MLstem(type[s], -1, -1, P);
471 new_fML = MIN2(energy, new_fML);
474 decomp = ggg[indx[j] + i] + n_seq * E_MLstem(0, -1, -1, P);
475 new_fML = MIN2(new_fML, decomp);
478 /* modular decomposition -------------------------------*/
479 for (decomp = INF, k = i+1+TURN; k <= j-2-TURN; k++)
480 decomp = MIN2(decomp, Fmi[k]+fML[indx[j]+k+1]);
482 DMLi[j] = decomp; /* store for use in ML decompositon */
483 new_fML = MIN2(new_fML,decomp);
485 /* coaxial stacking deleted */
487 fML[ij] = Fmi[j] = new_fML; /* substring energy */
491 int *FF; /* rotate the auxilliary arrays */
492 FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
493 FF = cc1; cc1=cc; cc=FF;
494 for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
497 /* calculate energies of 5' and 3' fragments */
501 case 0: for(j = TURN + 2; j <= length; j++){
503 if (c[indx[j]+1]<INF){
504 energy = c[indx[j]+1];
505 for(s = 0; s < n_seq; s++){
506 tt = pair[S[s][1]][S[s][j]];
508 energy += E_ExtLoop(tt, -1, -1, P);
510 f5[j] = MIN2(f5[j], energy);
514 if(ggg[indx[j]+1] < INF)
515 f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
518 for(i = j - TURN - 1; i > 1; i--){
519 if(c[indx[j]+i]<INF){
520 energy = f5[i-1] + c[indx[j]+i];
521 for(s = 0; s < n_seq; s++){
522 tt = pair[S[s][i]][S[s][j]];
524 energy += E_ExtLoop(tt, -1, -1, P);
526 f5[j] = MIN2(f5[j], energy);
530 if(ggg[indx[j]+i] < INF)
531 f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
537 default: for(j = TURN + 2; j <= length; j++){
539 if (c[indx[j]+1]<INF) {
540 energy = c[indx[j]+1];
541 for(s = 0; s < n_seq; s++){
542 tt = pair[S[s][1]][S[s][j]];
544 energy += E_ExtLoop(tt, -1, (j<length) ? S3[s][j] : -1, P);
546 f5[j] = MIN2(f5[j], energy);
550 if(ggg[indx[j]+1] < INF)
551 f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
554 for(i = j - TURN - 1; i > 1; i--){
555 if (c[indx[j]+i]<INF) {
556 energy = f5[i-1] + c[indx[j]+i];
557 for(s = 0; s < n_seq; s++){
558 tt = pair[S[s][i]][S[s][j]];
560 energy += E_ExtLoop(tt, S5[s][i], (j < length) ? S3[s][j] : -1, P);
562 f5[j] = MIN2(f5[j], energy);
566 if(ggg[indx[j]+i] < INF)
567 f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
578 #include "alicircfold.inc"
581 *** backtrack in the energy matrices to obtain a structure with MFE
583 PRIVATE void backtrack(const char **strings, int s) {
584 /*------------------------------------------------------------------
585 trace back through the "c", "f5" and "fML" arrays to get the
586 base pairing list. No search for equivalent structures is done.
587 This inverts the folding procedure, hence it's very fast.
588 ------------------------------------------------------------------*/
590 If s>0 then s items have been already pushed onto the sector stack */
591 int i, j, k, p, q, length, energy, up, c0, l1, minq, maxq;
596 length = strlen(strings[0]);
597 for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
598 type = (int *) space(n_seq*sizeof(int));
602 sector[s].j = length;
603 sector[s].ml = (backtrack_type=='M') ? 1 : ((backtrack_type=='C')?2:0);
606 int ss, ml, fij, fi, cij, traced, i1, j1, d3, d5, jj=0, gq=0;
607 int canonical = 1; /* (i,j) closes a canonical structure */
610 ml = sector[s--].ml; /* ml is a flag indicating if backtracking is to
611 occur in the fML- (1) or in the f-array (0) */
613 base_pair2[++b].i = i;
615 cov_en += pscore[indx[j]+i];
619 if (j < i+TURN+1) continue; /* no more pairs in this interval */
621 fij = (ml)? fML[indx[j]+i] : f5[j];
622 fi = (ml)?(fML[indx[j-1]+i]+n_seq*P->MLbase):f5[j-1];
624 if (fij == fi) { /* 3' end is unpaired */
631 if (ml == 0) { /* backtrack in f5 */
633 case 0: /* j or j-1 is paired. Find pairing partner */
634 for (i=j-TURN-1,traced=0; i>=1; i--) {
637 if (c[indx[j]+i]<INF) {
638 en = c[indx[j]+i] + f5[i-1];
639 for(ss = 0; ss < n_seq; ss++){
640 type[ss] = pair[S[ss][i]][S[ss][j]];
641 if (type[ss]==0) type[ss] = 7;
642 en += E_ExtLoop(type[ss], -1, -1, P);
644 if (fij == en) traced=j;
648 if(fij == f5[i-1] + ggg[indx[j]+i]){
649 /* found the decomposition */
650 traced = j; jj = i - 1; gq = 1;
658 default: /* j or j-1 is paired. Find pairing partner */
659 for (i=j-TURN-1,traced=0; i>=1; i--) {
662 if (c[indx[j]+i]<INF) {
663 en = c[indx[j]+i] + f5[i-1];
664 for(ss = 0; ss < n_seq; ss++){
665 type[ss] = pair[S[ss][i]][S[ss][j]];
666 if (type[ss]==0) type[ss] = 7;
667 en += E_ExtLoop(type[ss], (i>1) ? S5[ss][i]: -1, (j < length) ? S3[ss][j] : -1, P);
669 if (fij == en) traced=j;
673 if(fij == f5[i-1] + ggg[indx[j]+i]){
674 /* found the decomposition */
675 traced = j; jj = i - 1; gq = 1;
685 if (!traced) nrerror("backtrack failed in f5");
686 /* push back the remaining f5 portion */
691 /* trace back the base pair found */
694 if(with_gquad && gq){
695 /* goto backtrace of gquadruplex */
699 base_pair2[++b].i = i;
701 cov_en += pscore[indx[j]+i];
704 else { /* trace back in fML array */
705 if (fML[indx[j]+i+1]+n_seq*P->MLbase == fij) { /* 5' end is unpaired */
713 if(fij == ggg[indx[j]+i] + n_seq * E_MLstem(0, -1, -1, P)){
714 /* go to backtracing of quadruplex */
721 for(ss = 0; ss < n_seq; ss++){
722 tt = pair[S[ss][i]][S[ss][j]];
724 cij += E_MLstem(tt, S5[ss][i], S3[ss][j], P);
728 for(ss = 0; ss < n_seq; ss++){
729 tt = pair[S[ss][i]][S[ss][j]];
731 cij += E_MLstem(tt, -1, -1, P);
737 base_pair2[++b].i = i;
739 cov_en += pscore[indx[j]+i];
743 for (k = i+1+TURN; k <= j-2-TURN; k++)
744 if (fij == (fML[indx[k]+i]+fML[indx[j]+k+1]))
754 if (k>j-2-TURN) nrerror("backtrack failed in fML");
760 /*----- begin of "repeat:" -----*/
761 if (canonical) cij = c[indx[j]+i];
763 for (ss=0; ss<n_seq; ss++) {
764 type[ss] = pair[S[ss][i]][S[ss][j]];
765 if (type[ss]==0) type[ss] = 7;
769 if (cij == c[indx[j]+i]) {
770 /* (i.j) closes canonical structures, thus
771 (i+1.j-1) must be a pair */
772 for (ss=0; ss<n_seq; ss++) {
773 type_2 = pair[S[ss][j-1]][S[ss][i+1]]; /* j,i not i,j */
774 if (type_2==0) type_2 = 7;
775 cij -= P->stack[type[ss]][type_2];
777 cij += pscore[indx[j]+i];
778 base_pair2[++b].i = i+1;
779 base_pair2[b].j = j-1;
780 cov_en += pscore[indx[j-1]+i+1];
786 cij += pscore[indx[j]+i];
789 for (ss=0; ss<n_seq; ss++) {
790 if ((a2s[ss][j-1]-a2s[ss][i])<3) cc+=600;
791 else cc += E_Hairpin(a2s[ss][j-1]-a2s[ss][i], type[ss], S3[ss][i], S5[ss][j], Ss[ss]+a2s[ss][i-1], P);
793 if (cij == cc) /* found hairpin */
796 for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
797 minq = j-i+p-MAXLOOP-2;
798 if (minq<p+1+TURN) minq = p+1+TURN;
799 for (q = j-1; q >= minq; q--) {
801 if (c[indx[q]+p]>=INF) continue;
803 for (ss=energy=0; ss<n_seq; ss++) {
804 type_2 = pair[S[ss][q]][S[ss][p]]; /* q,p not p,q */
805 if (type_2==0) type_2 = 7;
806 energy += E_IntLoop(a2s[ss][p-1]-a2s[ss][i],a2s[ss][j-1]-a2s[ss][q],
808 S3[ss][i], S5[ss][j],
809 S5[ss][p], S3[ss][q], P);
812 traced = (cij == energy+c[indx[q]+p]);
814 base_pair2[++b].i = p;
816 cov_en += pscore[indx[q]+p];
823 /* end of repeat: --------------------------------------------------*/
825 /* (i.j) must close a multi-loop */
832 The case that is handled here actually resembles something like
833 an interior loop where the enclosing base pair is of regular
834 kind and the enclosed pair is not a canonical one but a g-quadruplex
835 that should then be decomposed further...
838 for(s=0;s<n_seq;s++){
842 mm += P->mismatchI[tt][S3[s][i]][S5[s][j]];
848 p < j - VRNA_GQUAD_MIN_BOX_SIZE;
850 if(S_cons[p] != 3) continue;
852 if(l1>MAXLOOP) break;
853 minq = j - i + p - MAXLOOP - 2;
854 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
855 minq = MAX2(c0, minq);
857 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
858 maxq = MIN2(c0, maxq);
859 for(q = minq; q < maxq; q++){
860 if(S_cons[q] != 3) continue;
861 c0 = mm + ggg[indx[q] + p] + n_seq * P->internal_loop[l1 + j - q - 1];
870 if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
871 minq = j - i + p - MAXLOOP - 2;
872 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
873 minq = MAX2(c0, minq);
875 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
876 maxq = MIN2(c0, maxq);
877 for(q = minq; q < maxq; q++){
878 if(S_cons[q] != 3) continue;
879 if(cij == mm + ggg[indx[q] + p] + n_seq * P->internal_loop[j - q - 1]){
888 for(p = i1 + 3; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
890 if(l1>MAXLOOP) break;
891 if(S_cons[p] != 3) continue;
892 if(cij == mm + ggg[indx[q] + p] + n_seq * P->internal_loop[l1]){
899 mm = n_seq*P->MLclosing;
901 for(ss = 0; ss < n_seq; ss++){
902 tt = rtype[type[ss]];
903 mm += E_MLstem(tt, S5[ss][j], S3[ss][i], P);
907 for(ss = 0; ss < n_seq; ss++){
908 tt = rtype[type[ss]];
909 mm += E_MLstem(tt, -1, -1, P);
912 sector[s+1].ml = sector[s+2].ml = 1;
914 for (k = i1+TURN+1; k < j1-TURN-1; k++){
915 if(cij == fML[indx[k]+i1] + fML[indx[j1]+k+1] + mm) break;
918 if (k<=j-3-TURN) { /* found the decomposition */
924 nrerror("backtracking failed in repeat");
927 continue; /* this is a workarround to not accidentally proceed in the following block */
931 now we do some fancy stuff to backtrace the stacksize and linker lengths
932 of the g-quadruplex that should reside within position i,j
935 int cnt1, cnt2, cnt3, cnt4, l[3], L, size;
938 for(L=0; L < VRNA_GQUAD_MIN_STACK_SIZE;L++){
939 if(S_cons[i+L] != 3) break;
940 if(S_cons[j-L] != 3) break;
943 if(L == VRNA_GQUAD_MIN_STACK_SIZE){
944 /* continue only if minimum stack size starting from i is possible */
945 for(; L<=VRNA_GQUAD_MAX_STACK_SIZE;L++){
946 if(S_cons[i+L-1] != 3) break; /* break if no more consecutive G's 5' */
947 if(S_cons[j-L+1] != 3) break; /* break if no more consecutive G'1 3' */
948 for( l[0] = VRNA_GQUAD_MIN_LINKER_LENGTH;
949 (l[0] <= VRNA_GQUAD_MAX_LINKER_LENGTH)
950 && (size - 4*L - 2*VRNA_GQUAD_MIN_LINKER_LENGTH - l[0] >= 0);
952 /* check whether we find the second stretch of consecutive G's */
953 for(cnt1 = 0; (cnt1 < L) && (S_cons[i+L+l[0]+cnt1] == 3); cnt1++);
954 if(cnt1 < L) continue;
955 for( l[1] = VRNA_GQUAD_MIN_LINKER_LENGTH;
956 (l[1] <= VRNA_GQUAD_MAX_LINKER_LENGTH)
957 && (size - 4*L - VRNA_GQUAD_MIN_LINKER_LENGTH - l[0] - l[1] >= 0);
959 /* check whether we find the third stretch of consectutive G's */
960 for(cnt1 = 0; (cnt1 < L) && (S_cons[i+2*L+l[0]+l[1]+cnt1] == 3); cnt1++);
961 if(cnt1 < L) continue;
964 the length of the third linker now depends on position j as well
965 as the other linker lengths... so we do not have to loop too much
967 l[2] = size - 4*L - l[0] - l[1];
968 if(l[2] < VRNA_GQUAD_MIN_LINKER_LENGTH) break;
969 if(l[2] > VRNA_GQUAD_MAX_LINKER_LENGTH) continue;
970 /* check for contribution */
971 if(ggg[indx[j]+i] == E_gquad_ali(i, L, l, (const short **)S, n_seq, P)){
973 /* fill the G's of the quadruplex into base_pair2 */
975 base_pair2[++b].i = i+a;
976 base_pair2[b].j = i+a;
977 base_pair2[++b].i = i+L+l[0]+a;
978 base_pair2[b].j = i+L+l[0]+a;
979 base_pair2[++b].i = i+L+l[0]+L+l[1]+a;
980 base_pair2[b].j = i+L+l[0]+L+l[1]+a;
981 base_pair2[++b].i = i+L+l[0]+L+l[1]+L+l[2]+a;
982 base_pair2[b].j = i+L+l[0]+L+l[1]+L+l[2]+a;
984 goto repeat_gquad_exit;
990 nrerror("backtracking failed in repeat_gquad");
997 /* fprintf(stderr, "covariance energy %6.2f\n", cov_en/100.); */
999 base_pair2[0].i = b; /* save the total number of base pairs */
1004 PUBLIC void encode_ali_sequence(const char *sequence, short *S, short *s5, short *s3, char *ss, unsigned short *as, int circular){
1007 l = strlen(sequence);
1011 /* make numerical encoding of sequence */
1012 for(i=1; i<=l; i++){
1014 ctemp=(short) encode_char(toupper(sequence[i-1]));
1019 /* use alignment sequences in all energy evaluations */
1024 ss[i] = sequence[i];
1027 ss[l] = sequence[l];
1044 if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.')) continue;
1048 for (i=1; i<=l; i++) {
1051 if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.')) continue;
1058 for(i=1,p=0; i<=l; i++){
1061 if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.'))
1064 ss[p++]=sequence[i-1]; /*start at 0!!*/
1069 for (i=l; i>=1; i--) {
1072 if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.'))
1080 PRIVATE void make_pscores(const short *const* S, const char **AS,
1081 int n_seq, const char *structure) {
1082 /* calculate co-variance bonus for each pair depending on */
1083 /* compensatory/consistent mutations and incompatible seqs */
1084 /* should be 0 for conserved pairs, >0 for good pairs */
1085 #define NONE -10000 /* score for forbidden pairs */
1088 int olddm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
1089 {0,0,2,2,1,2,2} /* CG */,
1090 {0,2,0,1,2,2,2} /* GC */,
1091 {0,2,1,0,2,1,2} /* GU */,
1092 {0,1,2,2,0,2,1} /* UG */,
1093 {0,2,2,1,2,0,2} /* AU */,
1094 {0,2,2,2,1,2,0} /* UA */};
1097 n=S[0][0]; /* length of seqs */
1099 if (RibosumFile !=NULL) dm=readribosum(RibosumFile);
1100 else dm=get_ribosum(AS,n_seq,n);
1102 else { /*use usual matrix*/
1103 dm=(float **)space(7*sizeof(float*));
1104 for (i=0; i<7;i++) {
1105 dm[i]=(float *)space(7*sizeof(float));
1107 dm[i][j] = (float) olddm[i][j];
1111 for (i=1; i<n; i++) {
1112 for (j=i+1; (j<i+TURN+1) && (j<=n); j++)
1113 pscore[indx[j]+i] = NONE;
1114 for (j=i+TURN+1; j<=n; j++) {
1115 int pfreq[8]={0,0,0,0,0,0,0,0};
1117 for (s=0; s<n_seq; s++) {
1119 if (S[s][i]==0 && S[s][j]==0) type = 7; /* gap-gap */
1121 if ((AS[s][i] == '~')||(AS[s][j] == '~')) type = 7;
1122 else type = pair[S[s][i]][S[s][j]];
1126 if (pfreq[0]*2+pfreq[7]>n_seq) { pscore[indx[j]+i] = NONE; continue;}
1127 for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
1128 for (l=k; l<=6; l++)
1129 score += pfreq[k]*pfreq[l]*dm[k][l];
1130 /* counter examples score -1, gap-gap scores -0.25 */
1131 pscore[indx[j]+i] = cv_fact *
1132 ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
1136 if (noLonelyPairs) /* remove unwanted pairs */
1137 for (k=1; k<n-TURN-1; k++)
1138 for (l=1; l<=2; l++) {
1139 int type,ntype=0,otype=0;
1141 type = pscore[indx[j]+i];
1142 while ((i>=1)&&(j<=n)) {
1143 if ((i>1)&&(j<n)) ntype = pscore[indx[j+1]+i-1];
1144 if ((otype<cv_fact*MINPSCORE)&&(ntype<cv_fact*MINPSCORE)) /* too many counterexamples */
1145 pscore[indx[j]+i] = NONE; /* i.j can only form isolated pairs */
1153 if (fold_constrained&&(structure!=NULL)) {
1154 int psij, hx, hx2, *stack, *stack2;
1155 stack = (int *) space(sizeof(int)*(n+1));
1156 stack2 = (int *) space(sizeof(int)*(n+1));
1158 for(hx=hx2=0, j=1; j<=n; j++) {
1159 switch (structure[j-1]) {
1160 case 'x': /* can't pair */
1161 for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
1162 for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
1170 case '<': /* pairs upstream */
1171 for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
1175 fprintf(stderr, "%s\n", structure);
1176 nrerror("unbalanced brackets in constraints");
1179 pscore[indx[j]+i]=NONE;
1183 fprintf(stderr, "%s\n", structure);
1184 nrerror("unbalanced brackets in constraints");
1187 psij = pscore[indx[j]+i]; /* store for later */
1188 for (k=j; k<=n; k++)
1189 for (l=i; l<=j; l++)
1190 pscore[indx[k]+l] = NONE;
1191 for (l=i; l<=j; l++)
1192 for (k=1; k<=i; k++)
1193 pscore[indx[l]+k] = NONE;
1194 for (k=i+1; k<j; k++)
1195 pscore[indx[k]+i] = pscore[indx[j]+k] = NONE;
1196 pscore[indx[j]+i] = (psij>0) ? psij : 0;
1198 case '>': /* pairs downstream */
1199 for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
1204 fprintf(stderr, "%s\n", structure);
1205 nrerror("unbalanced brackets in constraint string");
1207 free(stack); free(stack2);
1210 for (i=0; i<7;i++) {
1216 /*--------New scoring part-----------------------------------*/
1218 PUBLIC int get_mpi(char *Alseq[], int n_seq, int length, int *mini) {
1224 for(j=0; j<n_seq-1; j++)
1225 for(k=j+1; k<n_seq; k++) {
1227 for (i=1; i<=length; i++){
1228 if (Alseq[k][i]==Alseq[j][i]) ident++;
1231 if ((ident/length)<minimum) minimum=ident/(float)length;
1234 mini[0]=(int)(minimum*100.);
1235 if (pairnum>0) return (int) (sumident*100/pairnum);
1240 PUBLIC float **readribosum(char *name){
1248 int translator[7]={0,5,1,2,3,6,4};
1251 dm=(float **)space(7*sizeof(float*));
1252 for (i=0; i<7;i++) {
1253 dm[i]=(float *)space(7*sizeof(float));
1255 while(1) { /*bisma hoit fertisch san*/
1257 if (*line=='#') continue;
1259 i=sscanf(line,"%f %f %f %f %f %f",&a,&b,&c,&d,&e,&f);
1261 dm[translator[++who]][translator[1]]=a;
1262 dm[translator[who]][translator[2]]=b;
1263 dm[translator[who]][translator[3]]=c;
1264 dm[translator[who]][translator[4]]=d;
1265 dm[translator[who]][translator[5]]=e;
1266 dm[translator[who]][translator[6]]=f;
1275 PRIVATE void en_corr_of_loop_gquad(int i,
1277 const char **sequences,
1278 const char *structure,
1284 int pos, energy, en_covar, p, q, r, s, u, type, type2, gq_en[2];
1287 energy = en_covar = 0;
1289 while((pos = parse_gquad(structure + q-1, &L, l)) > 0){
1291 p = q - 4*L - l[0] - l[1] - l[2] + 1;
1293 /* we've found the first g-quadruplex at position [p,q] */
1294 E_gquad_ali_en(p, L, l, (const short **)S, n_seq, gq_en, P);
1296 en_covar += gq_en[1];
1297 /* check if it's enclosed in a base pair */
1298 if(loop_idx[p] == 0){ q++; continue; /* g-quad in exterior loop */}
1300 energy += E_MLstem(0, -1, -1, P) * n_seq;
1301 /* find its enclosing pair */
1302 int num_elem, num_g, elem_i, elem_j, up_mis;
1308 /* seek for first pairing base located 5' of the g-quad */
1309 for(r = p - 1; !pt[r] && (r >= i); r--);
1310 if(r < i) nrerror("this should not happen");
1312 if(r < pt[r]){ /* found the enclosing pair */
1319 /* seek for next pairing base 5' of r */
1320 for(; !pt[r] && (r >= i); r--);
1321 if(r < i) nrerror("so nich");
1322 if(r < pt[r]){ /* found the enclosing pair */
1325 /* hop over stems and unpaired nucleotides */
1326 while((r > pt[r]) && (r >= i)){
1327 if(pt[r]){ r = pt[r]; num_elem++;}
1330 if(r < i) nrerror("so nich");
1331 s = pt[r]; /* found the enclosing pair */
1334 /* now we have the enclosing pair (r,s) */
1337 /* we know everything about the 5' part of this loop so check the 3' part */
1339 if(structure[u-1] == '.') u++;
1340 else if (structure[u-1] == '+'){ /* found another gquad */
1341 pos = parse_gquad(structure + u - 1, &L, l);
1343 E_gquad_ali_en(u, L, l, (const short **)S, n_seq, gq_en, P);
1344 energy += gq_en[0] + E_MLstem(0, -1, -1, P) * n_seq;
1345 en_covar += gq_en[1];
1350 } else { /* we must have found a stem */
1351 if(!(u < pt[u])) nrerror("wtf!");
1355 en_corr_of_loop_gquad(u,
1364 en_covar += gq_en[1];
1368 if(u!=s) nrerror("what the ...");
1369 else{ /* we are done since we've found no other 3' structure element */
1371 /* g-quad was misinterpreted as hairpin closed by (r,s) */
1372 case 0: /*if(num_g == 1)
1373 if((p-r-1 == 0) || (s-q-1 == 0))
1374 nrerror("too few unpaired bases");
1379 for(cnt=0;cnt<n_seq;cnt++){
1380 type = pair[S[cnt][r]][S[cnt][s]];
1381 if(type == 0) type = 7;
1382 if ((a2s[cnt][s-1]-a2s[cnt][r])<3) ee+=600;
1383 else ee += E_Hairpin( a2s[cnt][s-1] - a2s[cnt][r],
1387 Ss[cnt] + a2s[cnt][r-1],
1392 for(cnt=0;cnt < n_seq; cnt++){
1393 type = pair[S[cnt][r]][S[cnt][s]];
1394 if(type == 0) type = 7;
1396 ee += P->mismatchI[type][S3[cnt][r]][S5[cnt][s]];
1398 ee += P->TerminalAU;
1402 energy += n_seq * P->internal_loop[s-r-1-up_mis];
1404 /* g-quad was misinterpreted as interior loop closed by (r,s) with enclosed pair (elem_i, elem_j) */
1408 for(cnt = 0; cnt<n_seq;cnt++){
1409 type = pair[S[cnt][r]][S[cnt][s]];
1410 if(type == 0) type = 7;
1411 type2 = pair[S[cnt][elem_j]][S[cnt][elem_i]];
1412 if(type2 == 0) type2 = 7;
1413 ee += E_IntLoop(a2s[cnt][elem_i-1] - a2s[cnt][r],
1414 a2s[cnt][s-1] - a2s[cnt][elem_j],
1425 for(cnt = 0; cnt < n_seq; cnt++){
1426 type = pair[S[cnt][s]][S[cnt][r]];
1427 if(type == 0) type = 7;
1428 ee += E_MLstem(type, S5[cnt][s], S3[cnt][r], P);
1429 type = pair[S[cnt][elem_i]][S[cnt][elem_j]];
1430 if(type == 0) type = 7;
1431 ee += E_MLstem(type, S5[cnt][elem_i], S3[cnt][elem_j], P);
1435 energy += (P->MLclosing + (elem_i-r-1+s-elem_j-1-up_mis) * P->MLbase) * n_seq;
1437 /* gquad was misinterpreted as unpaired nucleotides in a multiloop */
1438 default: energy -= (up_mis) * P->MLbase * n_seq;
1449 PUBLIC float energy_of_ali_gquad_structure( const char **sequences,
1450 const char *structure,
1459 short **tempS5; /*S5[s][i] holds next base 5' of i in sequence s*/
1460 short **tempS3; /*Sl[s][i] holds next base 3' of i in sequence s*/
1462 unsigned short **tempa2s;
1464 int *tempindx, *loop_idx;
1467 int en_struct[2], gge[2];
1469 if(sequences[0] != NULL){
1470 n = (unsigned int) strlen(sequences[0]);
1471 update_alifold_params();
1474 tempS = S; tempS3 = S3; tempS5 = S5; tempSs = Ss; tempa2s = a2s;
1475 tempindx = indx; temppscore = pscore;
1477 alloc_sequence_arrays(sequences, &S, &S5, &S3, &a2s, &Ss, 0);
1478 pscore = (int *) space(sizeof(int)*((n+1)*(n+2)/2));
1480 make_pscores((const short *const*)S, sequences, n_seq, NULL);
1483 pt = make_pair_table(structure);
1484 energy_of_alistruct_pt(sequences,pt, n_seq, &(en_struct[0]));
1485 loop_idx = make_loop_index_pt(pt);
1486 en_corr_of_loop_gquad(1, n, sequences, structure, pt, loop_idx, n_seq, gge);
1487 en_struct[0] += gge[0];
1488 en_struct[1] += gge[1];
1492 energy[0] = (float)en_struct[0]/(float)(100*n_seq);
1493 energy[1] = (float)en_struct[1]/(float)(100*n_seq);
1497 free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);
1499 /* restore old memory */
1500 S = tempS; S3 = tempS3; S5 = tempS5; Ss = tempSs; a2s = tempa2s;
1501 indx = tempindx; pscore = temppscore;
1503 else nrerror("energy_of_alistruct(): no sequences in alignment!");
1509 PUBLIC float energy_of_alistruct(const char **sequences, const char *structure, int n_seq, float *energy){
1515 short **tempS5; /*S5[s][i] holds next base 5' of i in sequence s*/
1516 short **tempS3; /*Sl[s][i] holds next base 3' of i in sequence s*/
1518 unsigned short **tempa2s;
1525 if(sequences[0] != NULL){
1526 n = (unsigned int) strlen(sequences[0]);
1527 update_alifold_params();
1530 tempS = S; tempS3 = S3; tempS5 = S5; tempSs = Ss; tempa2s = a2s;
1531 tempindx = indx; temppscore = pscore;
1533 alloc_sequence_arrays(sequences, &S, &S5, &S3, &a2s, &Ss, 0);
1534 pscore = (int *) space(sizeof(int)*((n+1)*(n+2)/2));
1536 make_pscores((const short *const*)S, sequences, n_seq, NULL);
1539 pt = make_pair_table(structure);
1540 energy_of_alistruct_pt(sequences,pt, n_seq, &(en_struct[0]));
1543 energy[0] = (float)en_struct[0]/(float)(100*n_seq);
1544 energy[1] = (float)en_struct[1]/(float)(100*n_seq);
1548 free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);
1550 /* restore old memory */
1551 S = tempS; S3 = tempS3; S5 = tempS5; Ss = tempSs; a2s = tempa2s;
1552 indx = tempindx; pscore = temppscore;
1554 else nrerror("energy_of_alistruct(): no sequences in alignment!");
1559 PRIVATE void energy_of_alistruct_pt(const char **sequences, short *pt, int n_seq, int *energy){
1563 energy[0] = backtrack_type=='M' ? ML_Energy_pt(0, n_seq, pt) : EL_Energy_pt(0, n_seq, pt);
1565 for (i=1; i<=length; i++) {
1566 if (pt[i]==0) continue;
1567 stack_energy_pt(i, sequences, pt, n_seq, energy);
1572 PRIVATE void stack_energy_pt(int i, const char **sequences, short *pt, int n_seq, int *energy)
1574 /* calculate energy of substructure enclosed by (i,j) */
1577 int *type = (int *) space(n_seq*sizeof(int));
1580 for (s=0; s<n_seq; s++) {
1581 type[s] = pair[S[s][i]][S[s][j]];
1587 while (p<q) { /* process all stacks and interior loops */
1591 if ((pt[q]!=(short)p)||(p>q)) break;
1593 for (s=0; s<n_seq; s++) {
1594 type_2 = pair[S[s][q]][S[s][p]];
1598 ee += E_IntLoop(a2s[s][p-1]-a2s[s][i], a2s[s][j-1]-a2s[s][q], type[s], type_2, S3[s][i], S5[s][j], S5[s][p], S3[s][q], P);
1601 energy[1] += pscore[indx[j]+i];
1603 for (s=0; s<n_seq; s++) {
1604 type[s] = pair[S[s][i]][S[s][j]];
1605 if (type[s]==0) type[s]=7;
1609 /* p,q don't pair must have found hairpin or multiloop */
1613 for (s=0; s< n_seq; s++) {
1614 if ((a2s[s][j-1]-a2s[s][i])<3) ee+=600;
1615 else ee += E_Hairpin(a2s[s][j-1]-a2s[s][i], type[s], S3[s][i], S5[s][j], Ss[s]+(a2s[s][i-1]), P);
1618 energy[1] += pscore[indx[j]+i];
1622 /* (i,j) is exterior pair of multiloop */
1623 energy[1] += pscore[indx[j]+i];
1625 /* add up the contributions of the substructures of the ML */
1626 stack_energy_pt(p, sequences, pt, n_seq, energy);
1628 /* search for next base pair in multiloop */
1631 energy[0] += ML_Energy_pt(i, n_seq, pt);
1637 PRIVATE int ML_Energy_pt(int i, int n_seq, short *pt){
1638 /* i is the 5'-base of the closing pair */
1640 int energy, tt, i1, j, p, q, u, s;
1649 do{ /* walk around the multi-loop */
1650 /* hop over unpaired positions */
1651 while (p < j && pt[p]==0) p++;
1653 /* get position of pairing partner */
1655 /* memorize number of unpaired positions? no, we just approximate here */
1658 for (s=0; s< n_seq; s++) {
1659 /* get type of base pair P->q */
1660 tt = pair[S[s][p]][S[s][q]];
1662 d5 = dangles && (a2s[s][p]>1) && (tt!=0) ? S5[s][p] : -1;
1663 d3 = dangles && (a2s[s][q]<a2s[s][S[0][0]]) ? S3[s][q] : -1;
1664 energy += E_MLstem(tt, d5, d3, P);
1671 energy += P->MLclosing * n_seq;
1673 for (s=0; s<n_seq; s++){
1674 tt = pair[S[s][j]][S[s][i]];
1676 energy += E_MLstem(tt, S5[s][j], S3[s][i], P);
1680 for (s=0; s<n_seq; s++){
1681 tt = pair[S[s][j]][S[s][i]];
1683 energy += E_MLstem(tt, -1, -1, P);
1688 energy += u * P->MLbase * n_seq;
1692 PRIVATE int EL_Energy_pt(int i, int n_seq, short *pt){
1693 int energy, tt, i1, j, p, q, s;
1701 do{ /* walk along the backbone */
1702 /* hop over unpaired positions */
1703 while (p < j && pt[p]==0) p++;
1704 if(p == j) break; /* no more stems */
1705 /* get position of pairing partner */
1707 for (s=0; s< n_seq; s++) {
1708 /* get type of base pair P->q */
1709 tt = pair[S[s][p]][S[s][q]];
1711 d5 = dangles && (a2s[s][p]>1) && (tt!=0) ? S5[s][p] : -1;
1712 d3 = dangles && (a2s[s][q]<a2s[s][S[0][0]]) ? S3[s][q] : -1;
1713 energy += E_ExtLoop(tt, d5, d3, P);